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SUMMARY 

This paper presents a computer program which has been developed for the flutter 
analysis, including the effects of rigid-body roll, pitch, and plunge, of swept-wing sub- 
sonic aircraft with a flexible fuselage and engines mounted on flexible pylons. The pro- 
gram utilizes a direct flutter solution in which the flutter determinant is derived by using 
finite differences, and the root locus branches of the determinant are searched for the 
lowest flutter speed. In addition, a preprocessing subroutine is included which evaluates 
the variable bending and twisting stiffness properties of the wing by using a laminated, 
balanced ply, filamentary composite plate theory. The program has been substantiated 
by comparisons with existing flutter solutions. 

The program has been applied to parameter studies which examine the effect of 
filament orientation upon the flutter behavior of wings belonging to the following three 
classes; wings having different angles of sweep, wings having different mass ratios, and 
wings having variable skin thicknesses. These studies demonstrated that the program 
can perform a complete parameter study in one computer run. The program is designed 
to detect abrupt changes in the lowest flutter speed and mode shape as the parameters 
are varied. 

INTRODUCTION 

The use of filamentary composite materials in aircraft structures offers a great 
potential for weight savings over conventional (all-metal) construction. Composites also 
introduce added versatility into the design process by allowing the structure to be better 
tailored to meet the design criteria. One such criterion is the prevention of flutter in a 
designated flight regime. The flutter speed can be quite sensitive to the distribution of 
structural stiffness. With the added versatility of composite material, the stiffness dis- 
tribution can be tailored to avoid flutter by controlling the filament orientation in each 
lamina of the structural wing box. 

Flutter is a complex phenomenon and the designer and analyst may have little 
intuitive feel for flutter even in all-metal airplanes since, in addition to the stiffness 



distribution, the flutter speed may be highly dependent upon the mass distribution, espe- 
cially the placement of engines, the flexibility of their supporting pylons, the rigid-body 
degrees of freedom of the plane, and the flexibility and mass of the fuselage. Further- 
more, it is the lowest flutter speed which is critical to the design and this speed can be 
significantly different in two seemingly similar designs since the lowest flutter speed 
may be associated with a different flutter mode in each design. Consequently, it wotild be 
desirable to have a computer program capable of analyzing the flutter behavior of 
composite -wing aircraft and determining the influence on the lowest flutter speed of com- 
posite structural parameters (such as filament orientation and lamina thickness) as well 
as other structural parameters (such as ermine mass and fuselage stiffness). 

The purpose of this paper is to present a flutter analysis program applicable to 
parametric studies of swept-wing aircraft and to apply it to study some of the flutter 
characteristics of composite wings. The program, which is written in FORTRAN IV, 
includes the effects of rigid-body roll, pitch, and plxmge and is applicable to both sym- 
metric and antisymmetric flutter modes. In addition, the effects of ermines mounted on 
flexible pylons and a flexible fuselage are incorporated into the program. 

The wing is modeled as a bending -twisting beam and the aerodynamic loads are 
modeled from modified swept-wing strip theory (refs. 1 to 3) in which finite aspect ratio 
and compressibility effects may be approximately accounted for. The program uses a 
direct flutter analysis approach in Which the flutter determinant is derived from energy 
considerations using finite differences. The root locus branches of the determinant are 
searched for the lowest flutter speed. 

The program contains a preprocessing subroutine which evaluates the variable 
bending and twisting stiffness properties of the wing by using a laminated, balanced ply, 
filamentary composite plate theory. The program also allows a parameter of interest 
(e.g., filament orientation, lamina thickness, and fuselage stiffness) to be automatically 
incremented over a user specified range, and the critical flutter speed is found at each 
parameter value, The program is demonstrated by studying the effect of filament orien- 
tation on the lowest flutter speed of wings belonging to the following three classes: 

(1) wings having different angles of sweep, (2) wings havii^ different mass ratios, and 
(3) wings having variable lamina thickness. 

The program has been verified by comparison with examples found in the literature 
(refs. 4 to 6) and with the SADSAM computer program (ref. 7). The difficulty of deter- 
mining the lowest flutter speed is illustrated by discussir^ an example reported in refer- 
ence 5. The present program predicted the same flutter speed as reference 5 but in 
addition found two lower flutter speeds. 
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SYMBOLS 


Aj. coefficients of superimposed displacement field as defined in equation (D3) 

^ll’^12’^22’^33 elements of elastic matrix as defined by equations (9) to (11) 

a, a local dimensional and nondimensional offsets of elastic axis from midchord, 

respectively (positive towards trailing edge), a = ^ 

b« 


^C’^C 


local dimensionai and nondimensional offsets of aerodynamic center from 

midchord, respectively (positive towards trailing edge), a^ = ^ 

°r 


reference value of semichord 


local dimensional and nondimensional semichords of airfoil normal to elastic 
axis, respeetively, b = 

b„ 

r 


local semichord of wing box 


Theodor sen* s circulation function 


local lift-curve slope 


complex flutter determinant 


dimensional and nondimensional fuselage diameter, respectively, d = - 

I 

Young's moduli parallel and transverse to filaments, respectively 


pylon and/or fuselage bendir^ stiffness 


reference value of wing bending stiffness 


EI,EI 


local dimensional and nondimensional wing bending stiffness, respectively. 


elements of flutter determinant defined by equation (D4) 
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composite shear modulus 


(GJ). 


pylon and/or fuselage torsional stiffnesses 


GJ,GJ 


local dimensional and nondimensional torsional stiffnesses of wii^, 


respectively, GJ = 


(ED, 


upward vertical displacement of any point along pylon 


h,h 


local dimensional and nondimensional upward vertical displacements of 
elastic axis, respectively, h = -^ 


Jacobian of flutter determinant defined by equation (D2) 


reduced frequency. 


bj.u) 


kn = 


cos A 


L,L 

I 


dimensional and nondimensional lengths of pylons, respectively, L = — 

br 

length of wing along elastic axis 


Mp bending moment on pylon at root 

reference mass per unit length along wing elastic axis 


m 


m,m 


me, me 


local dimensional and nondimensional mass distributions per unit length along 


elastic axis, respectively, m = 


m 

m- 


m^ 


dimensional and nondimensional mass of engines, respectively, m_ = 

mZ 


meJe,i»e‘^e dimensional and nondimensional rotational inertias of engine masses. 


respectively, m^J^ = 

nij.F 


N 


n\imber of finite -difference stations 
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N. 


number of composite reinforced laminas in each cover face 


^ll’^12’^22’^66 elements of reduced stiffness matrix defined in equations (12) 

to (15) 




local dimensional and nondimensional upward aerodynamic lifts per unit length 
of elastic axis on wings, respectively, Qi, = 

^ b^(EI)j. 


Q ,Q 


local dimensional and nondimensional aerodynamic torques on wings. 


respectively, Q = 


Q 


" (EI). 




dimensional and nondimensional upward vertical loads, respectively, applied 


to wings by pylons, qj^ = 


l>r(EI)r 


^m’*^m 


dimensional and nondimensional bending moments, respectively, applied to 

q 

wings by pylons, = -— 


q«,q« 


dimensional and nondimensional torques, respectively, applied to wings by 


pylons, q^ = 




“ (El), 


R 


ratio of wing length along elastic axis to reference semichord, — 




local dimensional and nondimensional radii of gyration, respectively, of wing 


r 

about elastic axis, 


kinetic energy of wings 


U 


torque applied to pylon at its root 
strain energy of wings 


u 


wing displacements along elastic axis 
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V 


airspeed 


w 

5W^,5We 

x,y,z 


a 

/3' 

e,e 

!j 

V 

0 

0 


K 


A 

X 


chordwise wing displacement, positive towards trailing edge 


upward vertical displacement of wing 


Virtual work of air loads and ermine masses, respectively 


aircraft coordinate system as shown in figure 1 


axial coordinate of pylon 


local dimensional and dimensionless offsets of wing center of gravity from 


X 


a 


elastic axis, x^ - ^ — 

’ a u, 


twist of wing, positive leading edge up 


quantity defined by equation (Ell) 


dimensional and nondimensional distances, respectively, between adjacent 

finite -difference stations, e = - 

I 

thickness coordinate of jth lamina 


coordinate of wing along elastic axis 


pitch about Y-axis 


ai^le of filament orientation, measured from ^ axis 

2 


reference wing mass ratio. 


rpb. 


m^ 


sweep ai^le 


nondimensional dynamic pressure. 


73t7-2 

pbj.i V 

(EI)r 
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ht’hl 

I 

p 

rj’ 7/1 

<t> 


Poisson’s ratios for the composite material 

chordwise coordinate of wing, positive towards trailing edge 

air density 

stress components in wing 
roll angle about X-axis 


0 ) 


nondimensional frequency, 



1/2 


circular frequency 


ANALYSIS 

Assumptions and Governing Equations 

Consider the model of a swept-wing subsonic aircraft with engines and flexible 
fuselage as shown in figure 1. The following assumptions concerning this model have 
been made; 

(1) All deformations are sufficiently small so that the governing dynamic equations 
are linear. 

(2) In-plane motions of the airplane are neglected; hence, only roll, pitch, and 
plunge of any point of the plane are allowed. 

(3) The tapered, variable -thickness wings are modeled as beams having an arbi- 
trary distribution of bending and twisting stiffnesses. 

(4) The engines are modeled as lumped masses (having rollii^ and plunging inertias) 
extended out in the plane of the wings mounted on flexible, massless pylons. The pylons 
are modeled as bending -twisting beams rigidly attached to the wings at the elastic axis. 

(5) A rigid carry -through beam joins the two wings and is rigidly attached to the 
elastic axes. 

(6) The fuselage is modeled as two flexible, massless, bending -twisting beams, one 
extending fore and the other aft of the wings. Each beam has a lumped mass at its tip 
with rolling and plunging inertia. The fuselage beams are rigidly attached to the carry- 
through beam. 
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The aerodynamic loads are provided from modified swept-wing strip theory as 
given in references 1 to 3 and repeated for completeness in appendix A. Air loads which 
do not act on the wings themselves are neglected. Due to the symmetry of the problem 
only symmetric and antisymmetric flutter modes need be Investigated; however, the 
limiting case of cantilever wing flutter is also considered. 

For harmonic motion, the lateral and torsional deflections of the wing may be 
expressed as 


h = 




( 1 ) 


a = 




( 2 ) 


where Sq and are complex amplitudes, w is circular frequency, and t is time. 
The governing equations which are valid for symmetric, antisymmetric, and cantilever 
wing harmonic motions are derived in finite-difference form in appendices B and C. They 
may be expressed as 


|[S] - - [A(A,n)]j[q] = 0 (3) 

where Qs^, and are the stiffness, mass, and aerodynamic load matrices, 
respectively. (The complex matrix does not include the mass effects of the air 
surroimdir^ the wing; these effects are incorporated into matrix (Iv^.) Also in equa- 
tion (3), 00 is the column of generalized displacements of the model, A is the dimen- 
sionless dynamic pressure, and n is the dimensionless complex frequency. 

For a given value of A, equation (3) represents a complex eigenvalue problem, for 
which a set of eigenvalues + iS2j and eigenvectors jqQ ®xist. Each pair of 

values (A,n) represent a root of the flutter determinant associated with equation (3) and 
the loci of all such roots constitute the root locus branches of equation (3) as shown in 
figure 2. Stable roots of equation (3) are those for which > 0. Consequently, flutter 
is indicated by a change in sign of from positive to negative and may be fotmd by 
tracing the root locus branches. 

Numerical Procedure for Tracii^ the Root Locus Branches 
and Determining Flutter 

The root locus branches are searched for flutter in the following manner; 

(1) At zero airspeed, A = 0, all the branches emanate from the natural frequencies 
of the system when the effective mass of the air is taken into accovint because 


8 



(4) 


lim[A(X,S2^ = 0 

For this reason the natural frequencies are found first; in the computer program, these 
frequencies are foxmd with a determinant plottii^ routine. 

(2) The branches are then traced by incrementing X from zero with the natural 
frequencies, determined in step (1), as starting points. 

(3) For each value of A, the complex eigenvalue O is found by driving the flutter 
determinant to zero by using a Newton-Raphson scheme whose iterative steps are given 
in appendix D, 

(4) A flutter instability is indicated by the vanishing of the imaginary part of the 
complex frequency; hence, each branch is traced out until flutter is indicated or a maxi- 
mum airspeed is attained. (As discussed in a subsequent section, traces need not nec- 
essarily start at A = 0.) 


Preprocessing the Wing Stiffness Distributions 

In this section, the bending and twisting stiffness distributions of the wing are pre- 
sented for a laminated, balanced ply, filamentary composite wing box as shown in 
figure 3. 

The cover of each box contains Nj^ layers and each layer contains an equal num- 
ber of filaments oriented in the plus and minus 0j directions (j = 1,2,. . .jN^). The 
filament orientation angle and the lamina thickness tj may be variable along the 
span and different in each lamina; however, complete symmetry is assumed to exist 
between the top and bottom covers of the box. 

It is assumed that the deformations of the wing obey the usual thin-plate Bernoulli- 
Euler assumptions and that chordwise bending may be neglected. On this basis, the 
bending and twisting stiffnesses of the wings, as derived in appendix E for cross sections 
nornial to the elastic axis, are 


El 



GJ 






a 


33,j 


(5) 

( 6 ) 
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where b* is the semichord of the structural wii^ box as shown in figure 3 and 
^o " 




( 8 ) 


in which Zrj, is the distance of the upper and lower covers of the box from the middle 
surface of the box, (See fig. 3.) In addition a^j^ j and a^g j are coefficients of the 
following elastic matrix which defines the orthotropic stress -strain law (ref, 8) in the 
jth lamina with filaments oriented in the plus and minus 0j directions: 



\ 



a, . 

— , 1 

jVf 

1 



'12, j 

0 


'12,3 

' 22,3 

0 


2a 


0 

0 

33,3 




From reference 8 the coefficients of the elastic matrix are given by 

^11 3 " ^ ^11 + ^(^12 ®3 ^3 ^22 ®3 

^33,3 = 2|(Qu + Q 22 - 2Qi 2 - 2Q66) 0. cos2 0. 

+ Q00^sin4 0. + cos4 0,jj 

where is the reduced stiffness matrix whose elements are 


(9) 


( 10 ) 


( 11 ) 


Qll = 




Q22 ~ 


E, 




Qi2 - 


^^t^t 
^ - "^t"t^ 


( 12 ) 


(13) 


(14) 
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^66 = 


(15) 


In equations (12) to (15), and are the Young’s moduli parallel and transverse to 
the filaments, respectively, and is the shear modulus. Equations (5) and (6) may be 
readily reduced to the case of an isotropic wing and to the case of a solid wing. 

COMPUTER PROGRAM 
Description 

A computer program denoted COMBOF (Composite Box Beam Flutter) has been 
developed for the flutter analysis of swept ^wing subsonic aircraft. The aircraft model 
contains a flexible fuselage and allows the engines to be mounted on flexible pylons. In 
addition the program takes into account the effects of rigid -body roll, pitch, and plunge. 

The user has the option of either supplying the bending and twisting stiffness dis- 
tributions of the wing or supplying the reduced stiffness matrix, filament orientation, and 
lamina thickness of a composite wing and allowing the program to automatically calculate 
the stiffness distribution of the wing by using a preprocessing subroutine governed by 
equations (5) to (15). When the preprocessor option is used, the program can perform a 
complete parameter study in one computer run, consecutively incrementing a parameter 
(e.g., fiber orientation, lamina thickness) over a range of interest, When calculating the 
flutter speed in such parameter studies, it may he possible to start from the flutter speed 
and frequency corresponding to the previous value of the parameter and to trace only the 
previous critical branch. Use of the previous lowest flutter speed in such a manner is 
reasonable so long as the lowest flutter speed or its associated frequency does not change 
discontinuously as the parameter is incremented. However, it may happen that as the 
parameter is incremented, a lower flutter speed may appear whose associated mode 
shape is different from that at the previous parameter value. This situation implies that 
the lowest flutter speed and/or its associated frequency have changed discontinuously 
with the parameter increment. In order to avoid such discontinuities, periodic searches 
of all the branches, starting from X = 0, are automatically performed by the program. 

The flow diagram of figure 4 shows the logic in COMBOF for finding the lowest 
flutter speed during a parameter study. Once the parameter of interest is chosen, the 
study proceeds as follows: 

(3) After the input of aircraft geometry, aerodynamics, and program control infor- 
mation, the program finds a user specified number of lowest natural frequencies. 

(j) A subroutine named RLF traces each branch, starting from the natural frequen- 
cies found in step(T), till a flutter speed is found or a specified maximum airspeed is 
obtained. 


11 



(J) The lowest flutter speed is chosen, or, if no flutter speeds have been found, the 
program may terminate. 

@ The parameter of interest is incremented, and the preprocessor computes new 
wing bending and twisting stiffness distributions. A count is kept so that after a specified 
number of increments a complete search of all branches is again initiated from X = 0. 
Consequently, discontinuities in the lowest flutter speed are cheeked for periodically 
during the parameter study. 

^ If the count has not been exceeded, subroutine RLF traces the branch on which 
the lowest flutter speed occurred for the previous parameter value while testing 

the first and second derivatives of X with respect to Oj along the 

branch to see if flutter is imminent 

the number of airspeed increments made along the branch 
the maximum airspeed limitation 

These tests may indicate that, for the parameter value, flutter is not likely to occur on 
the same branch as it occurred for the previous parameter value and a complete search 
of all branches is initiated. 

(6) The procedure continues until the specified number of parameter iterations is 
fulfilled or an internal decision is made to terminate the study. A description of the 
user input and options, an example computer output and a listing of the program are pro- 
vided in appendix F. 


Validation 

In order to verify the accuracy of the flutter analysis, several comparisons were 
made with examples found in the literature (refs. 4 to 6) and with the SADSAM (ref. 7) 
computer program. The results of these comparisons are given in tables I, n, and m. 

In order to agree with references 4 to 6, the aerodynamic coefficient 

dynamic center a^ were set equal to 2 it and -0.5, respectively. 

In table I, a imiform cantilever rectangular straight wing (with aspect ratio of 6.67) 
as given by Goland in reference 4 (with later correction of results in ref. 5) is considered. 
The table shows that COMBOF gives good agreement with Goland' s exact solution both in 
flutter speed and frequency even when as few as 10 finite difference stations are used. 

The SADSAM solution, although not doing as well in this example as COMBOF, also gives 
good results. The small differences between COMBOF or SADSAM and Goland's solution 
is probably due to the approximation of Theodor sen's circulation function used in both 
COMBOF and SADSAM instead of the exact Bessel Function form used by Goland. For all 
the calculations in table I, the branch on which flutter occurs is the first torsional branch; 


fc, and the aero- 

L “J 
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this is not meant to imply that only torsion is associated with the flutter mode. The 
designation of any root locus branch is derived from the designation of the natural fre- 
quency from which it emanates . 

An example of an aircraft having a rigid fuselage and in symmetric flutter is pre- 
sented in table II, The aircraft has viniform rectangular straight wings of aspect ratio 6.67 
and has weights attached to the wii^ tips. A solution to this problem is given by Goland 
and Luke in reference 5. They considered two cases in this investigation as described 
in the table. In each case, both COMBOF and SADSAM were in good agreement with ref- 
erence b for tile flutter speed and frequencies of the first torsional branch. However, in 
tracing out other root locus branches, COMBOF predicted that even lower flutter speeds 
existed for case 1 on the first and second bending branches, These speeds were later 
confirmed by SADSAM. This example serves to point out the ease with which the lowest 
flutter speed can be missed. 

Finally, a comparison with reference 6 for a uniform cantilever rectangular swept 
wing of aspect ratio 12.4 is presented in table III. Since reference 6 contains analytical 
and experimental results, comparisons were made with both. As in tables I and H, the 
comparisons show good agreement in all cases. 


DISCUSSION OF RESULTS FOR COMPOSITE -WING 
BOX PARAMETER STUDIES 

The program has been applied to parameter studies in which the effect of filament 
orientation upon flutter behavior has been examined for wings belonging to’ the following 
three classes; wii^s having different armies of sweep, wings having different mass 
ratios, and wings having variable lamina thicknesses. In each class, a cantilever wing 
was considered whose filament orientation angle was the same in each lamina and was 
uniform along the span (i.e., 0j = 0 = Constant), 

The studies were carried out by using the program logic shown in figure 4 
(previously described). For the cases presented herein the composite material was 
boron-epoxy and the appropriate moduli and other wing properties are given in table IV. 

In each case, the filament ai^le 9 was incremented in 101 steps from 0° to 90°. 
Among the relevant information output by the program was the variation of El and GJ 
with 6 . These data are shown in figure 5, Notice that while GJ is symmetric 
about 45°, El is asymmetric, decreasing with increasing 6 . 

The variation of the dimensionless dynamic pressure X with 9 is shown in fig- 
ure 6 for three swept wings. Values of X a.bove the flutter boundary are unstable, 
whereas those below are stable. The figure indicates that for 9 < 55° the effect of 
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sweeping the wing is to raise the flutter speed by a factor between l/ycosA and 
1/cos A. This effect is about the same as that recorded in references 6 and 9 for all- 
metal wings. However, when 9 > 55° sweeping the wing seems to have no definite effect 
on flutter. Moreover, in general, it appears that the effect of wing sweep is of secondary 
importance to the effect of fiber orientation, and, as anticipated, the maximum flutter 
speeds are attained when GJ is maximum or near 0 = 45°. This general shape of the 
curve bears much similarity to that presented in reference 10. 

As discussed in a previous section, COMBOF is designed to find the lowest value of 
the flutter parameter X (which corresponds to the lowest flutter speed) even if X is 
discontinuous. Such a situation occurs in figure 7 for a relatively high-mass -ratio wing 
(a heavy wing or low-air -density —high -altitude flight, m/mj,K = 64) as the filament 
orientation angle 9 is varied. The solid line in the figure represents the lowest flutter 
speed and in general follows the shape of GJ in figure 5. The curve contains two dis- 
continuities: one at 9 = 36° and one at 81°. However, the flutter boundary or bound- 
aries must be continuous and in the vicinity of such a discontinuity are usually multi- 
valued. COMBOF may also be used to predict these multivalued curves. (See appendix F 
for program option.) Hence the dashed curve in figure 7 completes the flutter boundary 
between 9 = 36° and 45°. 

Further details of the nature of the flutter boimdary in the vicinity of the disconti- 
nuity at 0 = 36° are shown in figure 8. The figure indicates the change in flutter mode 
from the first torsional branch (dashed curves) to the second bending branch (solid curves) 
as 9 increases from 35.1° to 45.0°. Flutter is indicated by a change in sign (from plus 
to minus) of the imaginary part of S2. Notice that at 35.1° the second bending branch 
almost goes unstable near X = 4 while at the same value of 9 the first torsional branch 
is unstable near X = 6.5. When 9 = 43.2° the second bending branch becomes unstable 
near X = 5, while the first torsional branch is also unstable but at X = 7.3, With further 
increase in 9 to 45°, the second bending branch is still unstable while the first torsional 
branch has become stable. 

Apparently the discontinuities in the curves of figure 7 are highly dependent upon the 
mass ratio of the wing. Figure 9 shows that, for mass ratios below that shown in figure 8, 
the discontinuity decreases as the mass ratio decreases. This trend is in agreement with 
results presented in references 11 and 12. When = 8 no discontinuity appears. - 

Again the curves generally follow the variation of GJ. 

The variation of X with 9 for variable thickness lamina is shown in figure 10. 

The thickness of the lamina varied linearly along the span from the root to the tip. The 
figure indicates that the discontinuity in X decreases in size as the lamina taper 
increases, and as anticipated, the flutter speed generally decreases with taper as well. 
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CONCLUDING REMARKS 


A computer program has been developed for the flutter analysis of composite swept^ 
wing subsonic aircraft with a flexible fuselage and engines mounted on flexible pylons 
including the effects of rigid-body roll, pitch, and plunge. The program uses a direct 
flutter solution in which the flutter determinant is derived by usii^ finite differences and 
the root locus branches of the determinant are searched for the lowest flutter speed. The 
air loads are generated internally by usii^ modified swept-wing strip theory, bi addition 
a preprocessing subroutine is included which can be used to evaluate the variable bendii^ 
and twisting stiffness properties of the wing using a laminated, balanced ply, filamentary 
composite plate theory. The program has been substantiated by comparison with existing 
flutter solutions. 

The program has been applied to parameter studies in which the effect of filament 
orientation upon flutter behavior has been examined for wings belonging to fhe following 
three classes: wii^s having different angles of sweep, wings having different mass 
ratios, and wings having variable lamina thicknesses. These studies demopstrated that 
the program can perform a complete parameter study in one computer run. Such studies 
make it possible to use the lolvest flutter speed, for one parameter value, in obtaining the 
lowest flutter speed at the next parameter value while making periodic cheeks for the 
occurrence of abrupt changes in the lowest flutter speed and/or frequency as the param- 
eter is varied. The occurrence of these abrupt chaises is associated with the lowest 
flutter speed shiftily from one flutter mode to another with a relatively small structural 
chaise. 

In addition, the studies Indicated that, as anticipated, the highest critical flutter 
speed for swept and unswept wings occurred when the filament orientation was near ±45°. 
However, for wings having high mass ratios, tiie critical flutter speed changed abruptly 
with small changes in filament angle when the filament orientation was near ±45°. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., May 22, 1974. 
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APPENDIX A 


AERODYNAMIC WING LOADS 

In this analysis the aerodynahiic loads are provided by modified strip theory 

(refs. 1 to 3) valid for constant amplitude harmonic motion of the form ei^t_ This 

theory represents an improvement of elementary strip theory in that it includes effects 

due to finite aspect ratio and compressibility. The theory employs a variable section 

lift-curve slope c, (rj) instead of 2jr, a variable section aerodynamic center aJrj) 

% 

instead of -1/2 (the quarter -chord), and a modified Theodorsen circulation function C(k) , 

which accomits for compressibility effects on the m^nitude of the lift and pitching 

moment. Appropriate spanwise distributions for c, ( 77 ) and a.(rj) are those for the 

c 

particular planform and Mach number of interest. (See refs. 1 to 3.) 

The lift and pitching moment per unit span at the ith finite-difference station may 
be expressed in dimensionless form as 






X 


n 


da. 

d?j 


tan A 


(Al) 
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APPENDIX A 



where the dimensionless (unbarred) loads are 

by 


^b,i " b' 


Q, . 1 ^ 

- h,r 


.(EI)j 




Q • = — ^ 
(El), 






(A3) 


in which T is the wing span (tip to root), bj. is a reference semichord and (El)j. is a 
reference bending stiffness. Also, 





and A is the sweep angle. Furthermore, 
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APPENDIX A 


- X cos^ 


A = pbj.r^v2 


cos^A 

(EI)r 


k =-^ 
“ cos A 

TTpb 2 

K = — 

mj. 


bj.w 

V cos A 


fi = a)r^^myy/(EI)j. 

An approximate form of Theodor sen’s circulation function, due to reference 7, is 
used in place of the exact Bessel function form; that is, 


^ , (W10.61ik„b,)(l^l.774uc„b,) 

^1 + 13.51iknbjj ^1 + 2.745iknb>^ 


(j = 1, . . N) (A4) 


Moreover, to correct for compressibility effects, may be modified in magnitude 

by the choice of an appropriate factor. (See ref. 3 for additional details.) 
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APPENDIX B 


GOVERNING EQUATIONS 


Since the aircraft is geometrically and structurally symmetric, only the symmetric 
and antisymmetric oscillations of the aircraft need be considered. Hence, the total 
elastic energy stored in both wings may be expressed as 



d?7 


(Bl) 


where a prime denotes differentiation with respect to rj, the coordinate aloi^ the elastic 
axis of the wing. The variation of eqiiation (Bl) may be expressed as 

N-1 

l6U=ie(Ei)iK^ah"H.€ 2 Mihl'aSi +|5(S)Nh>h'; 

i=2 

N 

+ ^ ^ (^‘^)i-l/2 “1-1/2 ®“i-l/2 
i=2 


where the elastic axis of the wing has been subdivided into N-1 equal segments of 
length 6 and the subscripts denote the finite -difference stations numbered as shown in 
figure 1. 

The kinetic energy of both wings may be expressed as 


where 



— - 2-2 
+ mrQ, a‘‘ 


At] 


(B3) 


m mass of wing per unit span 


r^ radius of gyration of cross section about elastic axis 

Xqj offset of center of gravity from elastic axis (positive towards the trailing edge) 
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and a dot signifies differentiation with respect to time. The variation of equation (B3) 
may be expressed as 




¥a,l 



- V«,N, 



(B4) 


The variation of the work done by the engines and fuselage may be expressed as 


N-1 

1 (%, A + V,i“i + 

i=2 

+ I ^(^h,N^^N ^m,N®^N ^o!,N®“n) 


where q, ,, q and q , are the vertical upward loads, torques, and moments 

{positive as shown in fig, 11) applied to the wings at the ith finite -difference station by the 
engines and fuselage. These loads are evaluated in appendix E, and are taken to be zero 
at the finite-difference stations where no engines are attached. 

The variation of the work done by the aerodynamic loads may be expressed as 
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N-1 




i=2 


|^(^h,N®^N '^Q!,N®®n) 


(B6) 


Substituting equations (B2), (B4), (B5), and (B6) into Hamilton's principle, 

>tr 


V ( 


5U - 5T - 6Wg - 6W^) dt = 0 


(B7) 


with harmonic motion of the form being assumed, and using the following finite - 

difference expressions for derivatives of Ej and oi. and their respective variations 


" ^{^i+1 ■ ^i-l) 


(B8) 




h’/ = 4(h.^i - 2h. + h,.,) 


(B9) 


“1-1/2 = |(“i - “i-l) 


(BIO) 


yield 


r^o fi 1 Q 1 (EI)^6hn r 

° 1 1_2 ■ ^^1 + ^o) ■*’4^ 'Im, J 13 Tr 2 + ^o) 


+ <EI)2(hj - 2h2 + hj) - l£4min2(hj - - |e^(Qh,i ^ %,i) 


1,4/ 


4'V,2|-3fcr - 21>1 *l>o) - 2(EI)2(h3 - 2h2 + hj) 


(EI)^6hj 




+ (EI)3(h4 - 2hg + h2) - 6 V2n2(h2 - - e4(Q^_2 + '‘h.2) 


(Equation continued on next page) 
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i— 3 


^ " ^^i+2 " + ^i) 


2(EDi(hi^l - 2hi + h._i) + (EI).^i(h. - 2h._j^ + h._2) - - o!jX^ j) 


“ ^ ^(^h,i ■*■ %,i) 2 ^ (^m,i+l ” ^m,i-l) 


(ED 6h. r 

Q®^N-2(%-l " -^-2 * ^N-s) 


- 2(EDj^_i(% ~ ^%-l %-2) + ■*■ ^N-l)’ 




- (vn-2 - 1 Vn)]-“^P^ •*■ 


‘N " “"N-1 ^N- 2 ) 


2h.T 1 + 


- (EOn^N+I - "n-i) - |'^”'n^^(''n - “n*«,n) - I'^Kk * «h,N) 


1,4/ 


.1, . _2h ^.K3_ 

2' Vn-^ e3fR2 *(2'“'n(“n+1 ^'^ + “n-i) J' Vn) e3fR2 


13 


[-«3fl3/2(“2 - “ 1 ) - - ^a,l\) ' |‘^('=>o/,l + «a,l) 




ef 


N-1 


i=2 




V. XA«.«Ui/ d . 

■^G%1/2 h+1 - «l) •»• <GJ>i-l/2 (“l - “i-l) - “^(“i^Q-,i ■ ^aA) 


(Eqviation continued on next page) 


22 



APPENDIX B 



In deriving equation (Bll) the variation of the kinetic energy has been integrated once by 
parts and the following dimensionless quantities have been Introduced: 



Not all the variations in equation (Bll) are independent of one another; some are 
constrained by the requirement of either symmetric or antisymmetric motion. 

Symmetric Modes 

For symmetric modes, 

<^)p = 0 (B12) 

where (p-^ the roll of the fuselage. However, since the wings are attached to the 
fuselage via a rigid carry-through beam, is also the roll of the wings at the first 

finite -difference station, hence 

0 - 4>-p = (p-^ = h* 2 ^ cos A + 0 ! sin A (B13) 


Expressing equation (B13) in finite-difference form, with dimensionless quantities intro- 
duced, gives 
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hQ = hg + 2eRaj^ tan A (B14) 

whose variation is 

6hQ = Shg + 2eRSa!j^ tan A (B15) 

Inasmuch as equation (Bll) must be valid for all independent variations of hj 
and 01^ which satisfy equation (B15), the governing finite-difference equations for sym- 
metric flutter may be obtained by setting the coefficients of each variation equal to zero. 
The resulting equations are 




0 = -(GJ)3/2(Q!2 “ “l) " I p ^ " hh) " ^a,l) 


+ 2 


(EI)j(h 2 - hj + eHOj tan a) + 


tan A 


eR 


(B16) 


0 = -2(EI)i(h2 + (EDg^hg - 2h2 + h^) - (h^ - 

■ 2 + %,l) + 2 ^^^m,2 “ 2€R(EI)^a!^ tan A 


(B17) 


e^mgO^ / 2 \ 

0 = -(GJ)5/2(Q!3 - 012 ) + (GJ)3/2(o! 2 " “l)" —3 — \^^a,2 “ %^2) 


~^%^a,2 + ^0!,2) 


(B18) 


C = 2(EI)j^^h2 - hj^ +€Rq!j^ tan A^ - 2(EI)2(hg - 2hg + hj^^ 

+ (EI)j(h^ - 2hg + hg^ - e^mgfi^^hg - a2^o!,2) 

-e^(Qh,2+'lh,3) ®1S) 
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0 = - «,) + (GJ)j.j/ 2(«, - «,.i) 

* »h,l) “ S.4- • -.N-2) 




0 = (EI)i^i(hj^2 - 2»1+1 * “i) - 2(EI)i(hi^i - 21>i + hi.ij 

+ (EI)i_i(hj - 2h,.j + hj.j) - e^m,n2(hi - 

- '^(«h,i 1h,i) - |'®(<lm,Ul - V,i+l) = 3,4,. . .,N-2) 


0 - ■(*^J)n-1/2(“n ■ “N-i) ^^'^)n-3/2 (“n-1 " “n-2) 

i^(“N-i^a,N-l " ^N-l^a,N-l) " ^^(^a,N-l ■*■ ^0!,N-l) 




0 - (EI)n-2(%-1 " ^^N-2 + ^^N-s) " “ ^%-l ■*■ *^N-2) 

~ 2 ^ (^m,N-2 " ^hIjN) 


0 = (GJ)n_i/2 {«n “ ®N-l) “ I ■ jj2 (“N^a,N “ *q;,N^n) 


(B20) 


(B21) 


(B22) 


(B23) 


(B24) 
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0 = " ^^N-1 ^N-2) ■ " “n^o!,n) 


2 ^^(^h,N ■*■ ^h,N) “ 2 ^^(V,N * ‘^mjN-l) 


- ~€^ I 


(B25) 


In addition the following equation, which resvilts from equation (Bll), was used to 
eliminate from equations (B23) ^d (B25) 


0 = (El)jj(hj^^j - 2hj^ + hj^_jj - 


(B26) 


Antisymmetric Modes 

In the antisymmetric case the constraints on the motion are 

0p = hp = 0 (B27) 

where ©p and hp are the pitch and plunge of the fuselage, respectively. Since the 
wings are rigidly attached to the fuselage j 

= (B28) 

where d is the fuselage diameter. Then as a consequence of equations (B27) and (B28) 
and assumption (6), the antisymmetric constraints take the following form* 


0 = 0p = Qfj cos A - (hg - h^j) 


(B28) 


hj = |<j|Raj sin A.+ (hj - h„) 


(B30) 


where 



I 

Equations (B29) and (B30) together imply 


(B31) 


(B32) 


26 



h 


1 ^(’'2 - \) 

1 4 € cos A 


APPENDIX B 


and their variation gives 


5 a 




tan A 
eR 


6 hj 


1 

4 e cos A 


(B33) 

(B34) 

(B35) 


Again, equation (Bll) must be valid for all independent variations of hj and a;| which 
satisfy equations (B34) and (B35); hence, the governing finite-difference equations for 
antisymmetric flutter are obtained as follows: 


0 = 


I (EI)]^ + i (EI)j I cos A + 1 (GJ)g /2 ^ “ | 


tan^ A 

"r^ 


+ l.m.x^,de2s22JaiU^^l 

16 ^ R cos A 8 


(EI)i /^v2 ^ j (El) 


1 


(El) 


2 d\ 


R cos A 8 cos2 aV } ^ cos A e 16 cos^ A 


2 2 ' 
1 ^ de SI , I ^ „ .2o2^ tan A 

m< — — +■ — m^x^, ^ e Si d — — 

32 ^ eos2 A 16 ^ R cos A 


ho + 


|Wi-j 


4 cos A £ 


- i (OJ), „ tan2 A + i ^ 1 

4 3/2 8 oiA 1 r 2 16 1 RcosA 

(El), ,, 1 (EDg d 1 mid2e2s22 


+ i ^^^^1 d _ 1 ^^^^2 /d\^ ^ 1 ^*^^2 d 


4 cos A ^ 8 cos^ A\^/ ^ ^ ^ 16 cos^ A ^ 62 cos^ A 


-^mid^x tan A 


^ C — 

16 ^ ^ C!,l cos A 




+ i e3/Q j + q + i ^ -/q i + q,^ o d cos A 

2 \ Q!jXy g COS A\ hjl iiyL I 0 1 X 1^2 


, tan A , 1 e d 


1 . 2 . 


. 1 (c* t\ ru tan A 
+ 2 ^2 -Jr 


(B36) 
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2 2 

0 = -(GJ)5/2(“3 - “2) ^^'^3/2 “2 " |{^2 " V ^|~(“2^ q !,2 " 


" ^^(^a!,2 + ^o!,2 ) 


(B37) 


0 = i (EDi^h^ . 1 f . h„J . 2(EI)2(h3 - 2h3 4 f ^ (EI)3(h, - 2h3 . h3) 


- €*A^(h2 - 02>!„,2) - * 2 '® ■‘m.S ' 


nirl 


- |((GJ)3/2[“2 - |(“2 - *H)) 




tan A _ 1 ,1 d\ ^2 ~ ^0 

eR 4 cos A 7 g2j^2 


,h -13!^*h3) 
I\ ^ 2 € cos A 


+ (EI)2(h„ - 2ho + j - i /l _d 

'‘^2^3 2 4 £ cos A / 4 1\2 6 co£ 




cos A 0!jl eR 


tan A 


(''2 - ''0) 




1 , 3 <- 


e cos A 


(B38) 


In addition, equations (B20) to (B25) are valid for antisymmetric modes. 

Cantilever Wii^ 

It is often assumed that the flutter behavior of an aircraft can be approximated by 
neglecting the rigid-body motion of the aircraft and any effect of fuselage flexibility. With 
these assumptions, the constraints on the fuselage are 

~ ~ ^ (B39) 
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By usii^ equations (B13), (B29), (B30), and (B39), the constraints take the form 

0!j = hj = 0 (B40) 

and 

h2 = hQ (B41) 

The wing is, therefore, cantilevered along a root normal to the elastic axis. Taking the 
variation of equations (B40) and (B41) yields 

6hj = 6hj = 0 (B42) 


and 


ehg = 6hQ (B43) 

Equations (Bll) must be valid for all independent variations of and hj^ which 
satisfy equations (B42) and (B43); hence, the governing equations for flutter of a cantilever 
wing are 


/ 2 \ 

0 = -(GJ)5/2(“3 - “2) + (0'’>3/2 “2 ’ "“2\“2''a,2 ‘ ='t<,2'»2) 

■ '^('^a ,2 + ‘10,2) 

0 = 2(EI)j hg - 2(EI)2(h3 - 2h2^ + (EI)3(h4 ^ 2h^ + hg) 

- e^n^mg^hg - XQj^gttg^ - + | ^^V,3 


In addition, equations (B20) to (B25) are valid for cantilever wing flutter. 
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EVALUATION OF ENGINE AND FUSELAGE LOADS 

A free -body diagram of the pylon -engine model is shown in figure 11. It is 
assumed that the pylon is rigidly attached to the wir^ at the ith finite -difference station. 
The bendii^ and twisting deformations of the pylon are assumed to be uncoupled and 
hence may be treated separately. 

Pylon Bending 

From elementary beam theory, bending deformations of the pylon are governed by 
d^h 

(EI)p = co^mgXphpCO) (Cl) 

where (EI)p is the bending stiffness of the pylon, Sg is the mass of die engine, and 
hp is the upward vertical displacement of the pylon, Integrating equation (Cl) twice 
yields 

hp = mgw2hp(0) ^ dx + ep(L) (l - Xp) + hpll) (C2) 

P 

where L is the length of the pylon. 

Evaluating equation (C2) at the ermine yields 


hp(0) = 


L6p(L) t hp(L) 


1 - nigW^ 


r^ x2 dx 
Jo (EI)p 


(C3) 


Then, from the equilibrium of the pylon. 


P = co^ni-. 
p e 




(C4) 
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and 


w^Lmg|^©p( l) - hp( L)j 

1 - 


(C5) 


where 

1 dx 

^^2" e Jo (EI)p 


(C6) 


and Pp is the downward vertical force at the pylon root and Mp is the moment at the 
pylon root. (See fig. 11.) For a uniform pylon, 



3(Epp 

rheC^ . 


(C7) 


Pylon Twist 

From elementary beam theory, torsional deformations of the pylon are governed 
by 



(C8) 


where nioJo represents the torsional inertia of the engine. Integrating equation (C8) 
and evaluating the result at the engine yield 



(C9) 


where 




dx 


- pL 

0 (GJ), 


(CIO) 


For a uniform pylon, 
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(GJ)p 

LnieJe 


(Cll) 


Hence from equilibrium, the torque at the pylon root is 


T„ = 




P' 


- (“yv ) 


(C12) 


Continuity of Displacements, Forces, and Moments at Pylon Root 

Enforcing the continuity of displacements and rotations at the ith finite -difference 
station where the pylon is rigidly attached to the wing’s elastic axis implies 

eJL) = COS A - sin A (CIS) 

<^p(l) = sin A + iij' cos A (Cl 4) 

hp(L) = -hj (C15) 


Moreover, from figvire 11 

’5h,i = 

eq^j j = Mp cos A + Tp sin A 
iqm i = Tp cos A - Mp sin A 


(C16) 

(C17) 

(C18) 


except for i = N where r/2 should replace r. Finally, substituting equations (CIS) to 
(CIS) into equations (C4), (C5), and (C12) and then substituting the resulting equation into 
equations (C16) to (Cl 8) yields, in dimensionless quantities, 






m. 




kcosA-$5%A)Ln.h, 
1 dr] R ' ^ 


(C19) 
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q Lf°SA . 

€r 2 \ ^ df7 R / 1 1 _ /n/j>„ 


Jl - 


^ J r 2 L sin a ^ ^ £2|i). _ »i ° - A . - ^ 
® ^ \ ^ *? / I - 


(C20) 


q . = J^ fi“R O', sin A + ^ 

'^m,i e e V i d?? R /i 


cos a\ cos a 


(n/nx)! 


r>2 („, „„„ A sin A W . u sin A 

eR dn R / jl _ fsi/S2^]2 


mg = 


2 

® ‘ (ECr 


„ 2 
"t 


In numerically evaluating the ei^ine loads, the following finite -difference equations were 
employed for symmetric flutter modes: 


-Roj tan A 


*"i+l - ^i-1 


(i=l) 


(i = 2,3,. . .,N-1) 


2 (El), 




X hj^ - (EI)j^hjj_j^ + i 


(i = N) 
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where 



(C23) 



(C24) 



(C25) 


and the subscript N designates the properties of an engine/pylon attached to the tip of 
the wing. The value of dhjq^dr; as expressed in equation (C22) was found from equa- 
tion (B26) and (C21). In addition 


dof. 


«2 " “l 


dn 




a- 1 - a. H 
1+1 1-1 


2e 


“n " “n-1 




For the case of antisymmetric flutter modes, 


d?j 


^2 “ ^0 
2e 


fh 1 ^(^2 ~ ^o) 

L ^ 4 e cos A - 


(i = 1) 


(i = 2,3,. . .,N-1) ) 


« = N,J 


(i=l)l 

( 1 = 2 ) 


(C26) 


(C27) 
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daj 

drj 




^ 


a = 1 ) 


(i = 2) 


(C28) 


J 


Equations (C22) and (C26) are valid for i > 2. For the case of a cantilever wing 


d?7 


^0 


1 

,2e 


(i=l) 


(i = 2) 


(C29) 




da 


d?7 






3 

2e 




(i=D 
(i = 2) 




(C30) 


For i > 2 the expressions for dh^yidtj and may be obtained from equations 

(C22) and (C26). 


Fuselage Loads 

Because of the similarity between the engine -pylon and fuselage models (compare 
assumptions (4) and (6)), the derivation of the fuselage loads may be readily obtained 
from equations (C19) to (C21) evaluated at i = 1. Consider the fuselage to be composed 
of two engine -pylon combinations, one extending forward and the other aft of the carry- 
throigh beam to which they are rigidly attached. Inasmuch as, for the case of symmetric 
or cantilever flutter modes, the roll, pitch, and plunge of the fuselage at the carry -throigh 
beam are the same as those at the first finite-difference station, equations (C19) to (C21) 
with i = 1 provide the fuselage loads. For antisymmetric flutter modes, the roll and 
pitch of the fuselage and the carry -throigh beam are still the same; however, the plvinge 
of the fuselage at the carry -through beam is zero, and that at the first finite -difference 
station is given by equation (B28). 
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EVALUATION OF THE FLUTTER DETERMINANT AND 
THE EXTRACTION OF ITS EIGENVALUES 

The value of the complex eigenvalue (frequency), U, at a value of X, is determined 
from a Newton-Raphson iteration techniqiie whose iterative steps are 


Hi ) Je,J 








!(3)i 





(Dl) 



8S2j dS2j 90jp^ 


(D2) 


and and (Djj j are the real and imaginary parts, respectively, of the flutter 

determinant Dj whose roots are sought on the jth iterative step. The partial deriva- 
tives appearii^ in equations (Dl) and (D2) are computed by finite -difference approxima- 
tions. To reduce the amount of computation a marching technique, subsequently 
described, is used to reduce the flutter determinant solved by equation (Dl) to order three. 


In this marching method, values of two displacements and one angle near the root of 
the wii^ are assumed. Through simultaneous application of all but three of the governir^ 
equations the displacements throughout the aircraft can be calculated from the assumed 
values. (The number of displacements and angles which must be assumed at the root is 
determined by the number of boundary conditions at the root and hence by the order of the 
governing equations, which is fourth order in h and second order in a.) The remaining 
three governii^ equations, valid at the wing tip, provide the reduced flutter determinant. 
The method proceeds as follows. 


Symmetric Flutter Modes 

For symmetric flutter modes, the following procedure is used: 

(1) Set hj^j^ = 1, h 2 ]^ = = 0, where the first subscript refers to the finite- 

difference station and the second to one of three assumed solutions 

(2) With the values set in step (1), solve for agi from equation (B16) 

(3) With the values set in step (1) and derived in step (2), solve for hgj^ from 

eqiiation (B17) 
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(4) Continue, using equations (B18) to (B28) to solve for hg^, . . 

^Nl’ ^Nl that order 

(5) Set t>22 “ ^ 

(6) Perform steps (2) to (4) calculatir^ a 22 > t^32> ®32’ ^42’ * * *’ ^N2’ “n 2 
(T) Set ^^3 ” ^23 " ^ and ^^3 “ ^ 

(8) Perform steps (2) to (4) calculatir^ agg, hgg, . . hj^g, a^jg 

The total solution is a linear superposition of the three assumed sets of displace- 
ments as 

3 3 

hi" ^ Aj.hj^j. oiy= ^ Aj.o!|j, (D3) 

r=l r=l 


The total solution must satisfy the three remaining governing equations (eqs. (B23) to 
(B25)). Hence at a value of = fij, one has 


fe<23) 

e(23) 

®2 

e(23^ 

®3 

p] 

e(24) 

®2 

e(24) 

®2 

e(24) 

®3 

Ag 

e(25) 

®3 

®2 

e(25) 

®3 

^3 

V-. 


where and e^^^^ (r = 1,2,3) are the complex values of the left-hand side 

of equations (B23), (B24), and (B25), respectively, for the rth set of trial solutions. For 
a nontrivial solution of equation (D4), the determinant of the coefficients must vanish; 
hence, the determinant of the matrix e is the reduced flutter determinant 

= l^ln=nj 


Antisymmetric Flutter Modes 

For antisymmetric flutter modes, the following procedure is used: 

(1) Set hgj = 1 and hgj^ = oigi “ ® 

(2) Solve for hQj^ from equation (B36), (Xg^^ from equation (B37), and h^j^ from 

equation (B38) 
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(3) Continue using equations (B20) to (B22) to solve for hgj, • • •, 

in that order 

(4) Set 1^22 ~ ^22 ~ ^ and ^32 “ ^ 

(5) Perform steps (2) to (4) 

(6) Set hgg = hgg = 0 and ags = ^ 

(7) Perform steps (2) to (4) 

The flutter determinant is then given by equation (D4). 

Cantilever Flutter Modes 

The same trial solutions are used for the cantilever flutter modes as for the 
symmetric case in conjunction with equations (B44), (B45), and (B20) to (B25). 


N1 


antl- 
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EVALUATION OF El AND GJ OF LAMINATED, BALANCED PLY, 
FILAMENTARY COMPOSITE WING BOX 

By employing the usual thin-plate Bernoulli-Euler assumptions and neglecting 
chordwise bendii^ of the wir^, the displacements w, u, and v take the form 




= h{n) - la{n) 


= -Z 


h'{^) - 



(El) 

(E2) 

(E3) 


Consequently, the strains are 

e„ = -zfh” - fo!’') 

^ dT] \ ' 


(E4) 


. = 1 (li + a = 

2 \a| 977/ 


(E5) 


€ 




T]Z 



(E6) 


If only the energy associated with the upper and lower covers of the box is considered, 
the strain energy per unit length along the elastic axis is 


where Zq is the distance from the middle surface of the box (z = 0) to die outer face of 
either the top or bottom cover of the box. Substituting equations (El) to (E3) into equa- 
tion (E4) yields 




dz d| 


(E8) 
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and substltutii^ equations (10), (11), (E4), (E5), and (E6) into equation (E8) yields. 


I ^ 2 "n 1 


i=i '•j-i 

Integrating through z gives 


'^r, = I.P S - 1“")' ^ 

3=1 


where 






and integrating through | gives. 


Nt 


= 1 ^j[2hl/('>”)^ ^ |(V)^o.")2 t 86*a33 j(o’)2] 


J=1 

Then introducii^ nondimensional quantities yields 

' 1^: 


= 


V 7 jj3 Aj ^ 

j=i L 


1 ^b* ^ + 4aoo 4b*R^ + |(b’^)‘'a, ^ . 

\d77 ) 3 


11,3^ 




‘33,J‘ 


*v3„ 1 6^ a} 




dn^jj 


(ElO) 


(Ell) 


(El 2) 


(E13) 


»b3|( 

where b has been replaced by bj.b . For high-aspect-ratio wii^s it is assumed that 
the last term on the right-hand side of equation (E13) may be neglected compared with 
the next to the last term (see ref . 13 for further explanation); hence. 


Nl ^ 


u 


1% 

3=1 

Furthermore, for a beam. 








(E14) 


U = -J— 

2rR3 


1_ / 2 \ 2 — , X • 

El / d^h\ ^i^2GJ/^Y 

[br \dr?2/ bj. \dr] ) 


(E15) 


Comparison of equations (Ell) and (E12) indicates that 
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El = 4b* 



a 


n,fy 


(E16) 


and 

- 

GJ - 16b 2 ^33,J^J 

3=1 


(E17) 


It is readily shown that equation (E16) is the same result derived from laminated plate 
theory and that equation (E17) reduces to Bredt' s formula for a thin-walled Isotropic tube 
with a rectai^ular cross section in which 


and 



- 1«1 

Zrp 



(E18) 


(E19) 
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COMPUTER PROGRAM 

The computer program COMBOF was written in FORTRAN IV on a SCOPE 3.1 sys- 
tem modified for Langley Research Center and executes and loads with a field lei^th 
specification of 60000 octal. The program is applicable to the flutter analysis (including 
rigid-body roll, pitch, and plunge) of composite swept -wing subsonic aircraft. To aid the 
user, a brief description of the input (which is supplied in Namelist format), ah example 
problem showing input and output, a program flow diagram (fig. 4), and a complete pro- 
gram listing are provided. 

Input Description 

NAMELIST /OPTION/ 

cantilever flutter modes 
symmetric flutter modes 
antisymmetric flutter modes 

option for tracii^ each root locus branch up to a maximum value of X 
(DAMMAX) while ignoring the occurrence of flutter 

normal running option for parametric studies in which the root locus 
branches are traced up to flutter or a maximum value of X 
(DAMMAX) 

flutter mode shapes are computed and printed 
flutter mode shapes are not computed 

maximum number of branches to be traced (also, number of natural 
frequencies calculated) 

one -dimensional array specifying order in which branches are to be 
searched (e.g., 2, 3, 1, 4, 5) where the branch number is deter- 
mined by the order of the natural frequency from which it emanates 

NAMELIST /TRACE/ 

DAM X increment used in tracing each branch 

IPROCD = 1 trace of branch will continue even if Newton-Raphson method fails to 
converge within specified number of iterations 


ISYM = 1 
= 2 
= 3 

KOPT s: 1 
= 2 

ISHAPE = 1 
= 2 

ILMAX 

IMORDER 
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IPROCD = 2 

ITMAX 

ISTOP 

DAMMAX 

DOMl 


search of branch terminates If Newton-Raphson method falls to con- 
verge within specified number of iterations 

maximum number of iterations for convergence of Newton-Raphson 
method when IPROCD = 2 

maximum number of X increments along any branch 

maximum value of X above which no search for flutter is made 

increment used in the finite difference evaluation of the partial 
derivatives appearing in the Newton-Raphson method 


NAMELIST /mOUT/ 


IX or JX = 1 
= 2 

IY= 1 

NAMELIST /PLAN/ 
SPAN = I 


display of two separately controlled sets of intermediate output 
no display of two separately controlled sets of intermediate output 
multiple parameter studies on one computer run 
single parameter study 


The following six symbols denote one -dimensional arrays which prescribe wing 
properties at each finite -difference station on the wing: 

A = a 

B = b 

SW ^ x^ 

RGSQ = rj 

BCH = 2B* 


ZT “ Zfji 

Q 


SWEEP 

N 

CL 


AC 


reduced stiffness matrix supplied by column 

sweep angle A of elastic axis in radians 

number of finite -difference stations along elastic axis 

one -dimensional array prescribing c, at each finite-difference 
station 

one -dimensional array prescribing a at each finite -difference 

station 
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NAMELIST /ENGINE/ 

The following six symbols denote one -dimensional arrays prescribing ei^ine, 
pylon, or fuselage properties of each engine or fuselage: 



ENGJ = 

ENJ = Je 
ENM = nig 
PYL = L 

log finite-difference stations at which pylons are attached. As discussed 

in appendix C, the fuselage is simulated by two engine-pylon com- 
binations, both having LOG = 1 but one extending fore and the other 
aft of the carry-through beam. 

NENG total number of ergines (including tixe number of fuselage lumped 

masses) 

NAMELIST /DESIGN/ 

EIREF = (EI)j. 

EMREF = mj, 

RHO = p 
RHOBAR 
BREF = bj. 

WTI 

LOW 


mass density of composite material 

one -dimensional array prescribing nonstructural weight per unit 
length of wing along elastic axis 

during a parameter study, periodic checks are made for discontinuities 
in lowest flutter speed by tracing all root locus branches starting 
from A, = 0; a count is kept of number of successive parameter 
increments for which no check is made and a check is initiated when 
this count equals integer value of LOW 


44 



APPENDIX F 


NL number of layers of composite material in top or bottom covers of 

wing box 

DB diameter of fuselage 

G acceleration due to gravity 

NAMELIST/ PARAM/ 

Two types of automatic parameter studies are available to the user. In the first 
option, ISTIPF = 1, the program internally calculates El, GJ, and m of the wing. 

The user supplies the engine properties in NAMELIST /ENGINE/ and the initial thick- 
ness distribution and filament orientation as functions of rj in a rectangular array AA; 
that is, 

t. = (AA)i^j 1 + t7(AA)2j + (AA)g^ . sin ^(AA)4^. (FI) 


= (AA)5^j 1 + 7?(AA)6,j 


+ (AA)(y j sin 



(F2) 


Equations (FI) and (F2) are employed in calculating El, GJ, and m. Ehiring the 
parameter study, specified elements of AA are incremented causii^ El, GJ, and m 
to be incremented. 


In the second option, ISTIFF = 2, the program increments engine, pylon, or fuse- 
lage properties and the user supplies and m in NAMELIST /STIFF/ and 

the initial engine-pylon or fuselage properties in array AA. (See listing of subroutine 
EVAL2 for details.) 

KVT number of elements of array AA which will be simultaneously 

incremented during study 

KV1,KV2 two one -dimensional arrays providing row and column designa- 

tion, respectively, of elements of AA which are to be incre- 
mented during study 


DEL one -dimensional array providing increment size of each element 

of AA specified by KVl, KV2 

MAX maximum number of parameter increments to be performed 

during study 

NAMELIST /STIFF/ (To be used only if ISTIFF = 2) 

EM,EI,GJ one -dimensional arrays which provide m, El, and GJ at each 

finite -difference station 
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NAMELIST /NATFQ/ 

OS starting value of used to search for natural frequencies of 

aircraft 

OE maximum value of during natural frequency search 

ODEL step size of used in natural frequency search 
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Program Listing 



PROGRAM COMBOFdNPUTf OUTPUT, TAPE5= INPUT, TAPE6=0UTPUT) 

A 

1 

c 


A 

2 

c 

COMPOSITE WING BOX FLUTTER 

A 

3 

c 


A 

4 

c 

NASA LANGLEY RESEARCH CENTER PROGRAM A3788 

A 

5 

c 

WING IN SUBSONIC FLIGHT 

A 

6 

c 


A 

7 

c 

RHO^AIR OERSITY • 

A 

8 

c 

WING PLANFGRM - INPUT 

A 

9 

c 


A 


e 

SPAN 

A 

11 

c 

A=OFFSET OF ELASTIC AXIS FROM SEMI-CHdRD 

A 

12 

c 

BOX SEAM CHORD - BCH(RODAL POINT) 

A 

13 

c 

RGSQ^'RACIUS OF GYRATION ABOUT ELASTIC AXIS 

A 

14 

c 

SW=OFFSET OF C.G. FRCM ELASTIC AXIS 

A 

15 

c 

EM=WING MASS PER UNit SPAN 

A 

16 

c 

EI=W1NG BENDING STIFFNESS 

A 

17 

c 

GJ=WING TORSIONAL STIFFNESS 

A 

18 

c 

EIR£F,EMREF=R£FeRENCE BENDING STIFFNESS AND WING MASS RESPECTIVELY 

A 

19 

c 

BREF*R£FERENCE WING SEMI-CHORD 

A 

20 

c 

ZT= ONE half the DISTANCE BETWEEN TOP ANO BOTTOM COVERS 

A 

21 

c 

RHOBAR=WING BOX COMPGSITE MATERIAL DENSITY 

A 

22 

c 

NL=NUMBER OF LAMINAS IN TOP OR BOTTOM WING BOX COVERS 

A 

23 

c 

WTI=W£IGHT OF WING PER UNIT SPAN INOT INCLUDING COMPOSITE" WEIGHT! 

A 

24 

c 

G*ACCELERATION DUE TO GRAVITY IN APPROPRIATE UNITS 

A 

25 

c 

D8-OIAMETER OF FUSELAGE 

A 

26 

c 

Q^REDUCED STIFFNESS MATRIX 

A 

27 

c 


A 

28 

c 

ENGINE DATA - INPUT 

A 

29 

c 

ENEI=B6NDING STIFFNESSES OF PYLONS 

A 

30 

c 

ENGJ=TORSIONAL STIFFNESSES OF PYLONS 

A 

31 

c 

ENJ=TORSIONAL INERTIA OF ENGINE 

A 

32 

c 

NENG=NUM8ER OF ENGINES (INCLUDE FUSELAGE IN COUNT) 

A 

33 

c 

LOC=NGOAL LOCATION OF ENGINE 

A 

34 

c 

PYL=PYLON LENGTH 

A 

35 

c 

ENM=ENGINE MASSES 

A 

36 

c 


A 

37 

c 

MATRIX AA CONTAINS INFORMATION RELATING TO THE LAMINA THICKNESS 

A 

38 

c 

AND FILAMENT ORIENTATION WHEN ISTIFF=1 AND RELATES TO ENGINE/PYLON 

A 

39 

c 

INFORMATION WHEN ISTIFF=2. ( SEE EVALl AND EVAL2) 

A 

40 

c 


A 

41 


COMMON // IXfJXt I£RRCR,KOPT/CUM/OM,LAMBDA,RHO,EM160) ,E I( 61 ) ,G J ( 61 ) 

A 

42 


1,A(60),8(60),SWI60),AR,RGSO(6D),MU,E(3,3) ,W(60) ,AL(60> ,NN, SWEEP, IP 

A 

43 


2ROCO,ITMAX,ISYM,CL(90I ,AC(90)/STATIC/IST, ALO,CMAC, EL, ISHAPE/PYLON/ 

A 

44 


3BPY(lO),CPYaO),OPY( i0J,ENJ(i0),ENM{10),ENEia0) ,ENGJ( 10) ,PYL(10) , 

A 

45 


4L0CU0) ,NENG/BLK1/RH0BAR,SPAN,BCH(60) ,G/BLK2/0 ( 3,3 ) , ZT (60 ) , EIREF, E 

A 

46 


5MREF,0BEGIN,LBEGIN/FREQ/0M1( 5),DOMl,OAM,OMMAXl,toAM,OAMMAX 

A ^ 

47 


DIMENSION AA(8,4), IM0RDER(5), WTK60), XVK4), KV2(4), DEL(4) 

A 

48 


REAL LAMBDA, LBEGIN 

A 

49 


COMPLEX OM,W,AL,E,OBEGIN 

A 

50 
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c 


A 

51 

c 

INOUT DATA IS IN NAMELIST FORM 

A 

52 

c 


A 

53 


NAMELIST /OPTION/ I SYM,KGPTf I SHAPE » ILMAX f IMORDER 

A 

54 


NAMELIST /TRACE/ DAM, IPROCO.ITMAX, ISTOP,OAMMAX,DOMl 

A 

5 5 


NAMELIST /INOUT/ [X,JX,IY 

A 

56 


NAME LI S T /PLAN/ S PAN , A , B , SW , R G SQ , B C H , Z T » Q , S WE E P , N » CL , AC 

A 

57 


NAMELIST /ENGINE/ ENEI ,ENGJ,ENJ,ENK,PYL,LOC,NENG 

A 

58 


NAMELIST /DESIGN/ EIREF,EMREF,RH0,RH0BAR,BREF,WTI,L0W,NL,DB,S 

A 

59 


NAMELIST /PARAM/ I STIFF , AA,KVT, KVl ,KV2 ,OEL,MAX 

A 

60 


NAMELIST /STIFF/ EM,EI,GJ 

A 

61 


NAMELIST /START/. OBEGINtLBEGIN 

A 

62 


NAMELIST /NATFQ/OSfOEfOOEL 

A 

63 

1 

READ OPTION 

A 

64 


READ TRACE 

A 

6 5 


INOUT “ 

A 

66 


READ PLAN 

A 

67 


READ ENGINE 

A 

68 


READ DESIGN 

A 

69 


READ PARAM 

A 

70 


IF (ISTIFF.EQ.2I READ STIFF 

A 

71 


( iStiFF. EQ.1. AND. NL.EQ.O) READ STIFF ~ ~ 

A 

72 

c 


A 

73 

c 

FOR FLUTTER SET iFORC==l 

A 

74 


I FORD =1 

A 

75 


IST=IFORD 

A 

76 


IM=0 

A 

77 


rSTART=3 

A 

78 


DBAR=.5*OB 

A 

79 


MU=N 

A 

80 


EL=N-1 

A 

81 


L8EGIN=0. 

A 

82 


0BEGIN=0* 

A 

83 


IF ( ISTART.EQ.2) READ START 

A 

84 


READ NATFQ 

A 

85 


PRINT 26 

A 

86 


PRINT 34 

A 

87 


PRINT 27, IX, JX,N, I FORD, MAX, SPAN, SWEEPS (KVK t ) .KV2C I ) , DEL ID , I== 1, K 

A 

88 


IVT \ 

A 

89 


PRINT 41, IPROCD,ITMAX,ISTOP,lSHAPi,ISYM, ISTART, ILMAX. KOPT 

A 

90 


PRINT 38 

A 

9l_ 


PRINT 28, (I,Am ,8m,SW(I),RG$Qin,I = l,MU) 

A 

92 


PRINT 49 

A 

93 


PRINT 50, (I ,BCH( n ,ZT( I),WTI(I ),I=1,MU) 

A 

94 


PRINT 29, {(Q(I ,J) , ,I-1,3» 

A 

95 


PRINT 30, COMI,DAM,DAMMAX 

A 

96 


PRINT 32, RH0tRH08AR,EIREf,EMREF,BREF 

A 

97 


IF (NENG.Ee.OI GOTO 2 

A 

98 


PRINT 31, (LOCI I) ,ENEI( U ,ENGJ(H ,ENJ(I » ,ENM(i ) , PY L( I) , 1 = 1 , NENG • 

A 

99 


GO TO 3 

A 

100 

z 

PRINT 43 

A 

101 

3 

IF nSTIFF.EQ,2) GO TO 5 

A 

102 


IF (NL,EO.O) PRINT 47 

A 

103 


IF (NL.EQ.OJ GO TO 7 

A 

104 


PRINT 44 

A 

105 


DO 4 I =1 , 8 

A 

106 

4 

PRINT 33, (AA(I, J),J=l,NLy 

A 

107 


GO TO 7 

A 

108 

5 

PRINT 45 

A 

109 


00 6 1=^1, 5 

A 

110 

6 

PRINT 33, ( AAH, Jl , J = 1,NENG) 

A 

111 

7 

CONTINUE 

A 

112 


IF ( ISTART.EQ.2I PRINT 40, OBEGIN,LBEGIN 

A 

113 


IERR0R=0 

A 

114 


IERRS-0 

A 

115 
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c 


A 

116 

c 

NON-Dl'MENSiONAUI Ifc 

A 

117 

c 


.. ^ 

118 


RBI?EE = l./BREf " ' “ 

A 

119 


AR=SPAN*RBREF 

A 

120 


REIREF=1./EIREP 

A 

121 


UBAR^OBAR/SPAN 

A 

122 


RHO=3.141592658979*BP.EF*'6RiEF*RHd/EMREF 

A 

123 


REMREF=1. /EMREF 

A 

124 


DO 8 1=1, MU 

A 

125 


BII)=B(I)*RBREF 

A 

126 


A(I) = Am*RBREF 

A 

127 


SW(I)=SW<I)*RBRE”F 

A 

128 

8 

RGSQ( n=RGSQ(I )*RBREF*RBREF 

A 

129 


PRINT 35 

A 

130 


PRINT 38 

A 

131 


PRINT 28 , (I,A(n,B(I),Swm,RGSQm,I«l,MU) 

A 

132 


IF (NENG.EQ.O) GO TO 10 

A 

133 


DO 9 1=1 ,NENG 

A 

134 


ENEK II=ENEim*REIREF 

A 

135 


ENGJM )=ENGJ(n^REIR£F 

A 

136 


ENJ( n=ENJ(I »/(SPAN*SPAN) 

A 

137 


ENM( n=ENMn )*REMR£F/SPAN 

A 

138 

9 

PYL( n=PYLm*RBREF 

A 

139 _ 


PRINT 31 , (LOG ( 1 1 , ENEI ( 1 1 , E NG J ( 1 1 , EN J ( I) , ENM C II , PYL ( H , 1^1 , NENG ) 

A 

140 

10 

IF ( ISTIFF.EQ.I.AND.NL.GT.OJ GO TO 12 

A 

141 


DO 11 1=1, MU 

A 

142 


tl ( I ) =EI( I )<^R£IREF 

A 

143 


GJ { I ) = GJ( n^'REIREF 

A 

144 

11 

EMm=EMm«REMREF 

A 

145 

12 

CONTINUE 

A 

146 


KKAPPA=1* ARHO 

A 

147 


PRINT 39, ARfRKAPPA 

A 

148 


PRINT 51 

A 

149 


AR=AR/EL 

A 

150 


D0M1=D0M1/EL^<'2 

A 

151 


0AM=DAM/EL**3 

A 

152 


OBAR=OBAR*EL 

A 

153 

C 

PRE CALCULATIONS FOR ENGINES 

A 

154 


IF (NENG.EQ.O! GO TO lA 

A 

155 


DO 13 IJ=1,NENG 

A 

156 


ENM( IJ)=ENM(IJ)*EL 

A 

157 


ENJ(IJI=ENJ(IJ)*Et*EL 

A 

158 

C 

l/BPY = FIRST BENDING FREQ. OF PYLON 

A 

159 

C 

1/DPY = FIRST TORSIONAL FREQ. OF PYLON 

A 

160 


BPY( lJI = PYLnj)*PYL( U l*PYL( U )*ENM( IJ )/ ( 3.*ENEI { IJ 1 »/(AR*AR*AR) 

A 

161 


CPY(tJ)=ENJ( IJ>*ENM( IJ» ~ 

A 

162 

13 

DPY( t JJ=CPY( IJ)*PYL( IJ»/ (ENGJ( I J)*AR) 

A 

163 

14 

CONTINUE 

A 

164 

C 


A 

165 

C 

CALL EVALl TO PREPROCESS EI,GJ AND EM IF IST!FF=1 

A 

166 

C 

CALL EVAL2 TO PREPROCESS ENGINE DATA IF ISTIFF=2 

A 

167 

C 


A 

168 


IF { ISTIFF.EQ.I! CALL EVALl IAA,WTI,NU ~~ 

A 

169 


IF nSTlFF,EQ.2) CALL EVAL2 (AA» ~ 

A 

170 


PRINT 36 “ 

A 

171 
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MU1=MU«-1 

A 

172 


PRINT 37, a,Eim,Gjin,EMm,cLm,Acm,i=i,Mui 

A 

173 


IF IISfART.NE.3) GO TO 15 

A 

174 


CALL NATRL ( I LHAX ,OS ,OE,ODEL 1 

A 

175 


PRINT 42, (OMl(n,l=l, ILMAX) 

A 

176 


IM=0 

A 

177^ 

15 

CALL RLF tOBEGIN,LBEGIN, ISTART,ILMAX,IM,IMOROER) 

A 

178 


V=S0RT(LBEGiN*EIREF/ISPAN**3)/BREF*3.141592658979*(BREF«BReF»/{EMR 

A 

179 


lEF*RHOn 

A 

180 


IF ( lERROR.EQ.OI PRINT 46, V 

A 

181 


PRINT 51 

A 

182 


IT=1 

A 

183 


KX = 1 

A 

184 


If IM AX. EQ.O) GO TO 25 

A 

18 5 

C 


A 

186 

C 

INCREMENT PARAMETER 

A 

187 

c 


A 

188 


DO 24 KZ-lf MAX 

A 

189 


IF (lERRS.GT. 2. AND. KOPT.EQ.il GO TO 24 

A 

190 


IERRS=0 

A 

191 


CALL VARY ( AA.KVl ,KV2, KVT,DEL,wtl ,NL , I STIFF 1 

A 

192 


IF ( ISTIFF,eQ.2l PRINT 48, LLOC Til , ENE I III , ENGJ ( 1 1 , EN J ( I 1 , ENMI I) ,P 

A 

193 


1YUII,I=1,NENGI 

A 

194 


If nSTIFF.EQ.l) PRINT 36 

A 

195 


IF ilSTIFF.EQ.il PRINT 37, M , El 1 1 1 ,G JM 1 ,EM(I 1 ,CL ( M , AC 1 1 1 , 1=1 ,MU 

A 

196 


1) 

A 

197 


If (ISTIFF.E0.2) GO TO 17 

A 

198 


PRIfSiT 44 

A 

199 


00 16 1=1,8 

, A 

200 

16 

PRINT 33, I AA( I ,J) , J=1,NLJ 

A 

201 


GO TO 19 

A 

202 

IT 

PRINT 45 

A 

203 


DO 1 8 1=1 ,5 

A 

204 

18 

PRINT 33, ( AA( I, J I , J=l,NeNGI 

A 

205 

IS 

CONTINUE 

A 

206 


IF (KX.LT.LOW) GO TO 20 

A 

207 


ISTARt=3 

A 

208 


KX=0 

A 

209 

2C 

IF (I START. £0.21 GO TO 22 

A 

210 

21 

XALL NATRL 1 ILMAX ,0S , OE, ODEL 1 

A 

211 


PRINT 42, lOMlII 1 ,I=i,ILMAXI 

A 

212 

22 

CALL RLF (Q8EGIN,LBEGIN, ISTART, ILMAX, IM,IM0RDER1 

A 

213 


V=SQR TIL 8EGIN*EIREF/|SPAN**3I/BREF*3. 141592658979* (BREF*BREFI/|EMR 

A 

214 


1EF*RH01I 

A 

215 


PRINT 46, V 

A 

216 


PRiNt 51 

A 

217 


KX=KX+1 

A 

218 


If ( IOPT.EQ.il GO TO 23 

A 

219 

23 

IERRS=IERRS+IERROR 

A 

220 


IF IIERRS.GT.2) GO TO 24 

A 

221 


IF ( lERROR^GT.O) GO TO 21 

A 

222 

24 

CONTINUE 

A 

22 3 

25 

IF (lY.EQ.l) GO TO 1 

A 

224 


STOP 

A 

22 5 _ 
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c 


A 

226 

26 

FORMAT aHl,35Xf6HCOMBOFf3X,28H(COMPOSITE WING BOX FLUTTER I , /, 35Xf 

A 

227 


l42(lH*lt///) 

A 

228 

27 

FORMAT (10X,3HIX=f 11 1 5X ,3HJX= , 1 1 , 5X, 2HN= , 1 2, 5X f 6HI F0RD= , I2,5X,4HMA 

A 

229 


1X=, I3,5X,5HSPAN='. E16.8,6X,6HSWEEP = ,E16.8//t <10X,llHKVltKV2,0EL, 2X, 

A 

230 


2I2»3X,I2,3X,E16.8),/) 

A 

231 

23 

FORMAT (5X,I2,6X,4(E1A.6,4X)) 

A 

232 

29 

FORMAT (//,15X,24HRECUCEC STIFFNESS MATR IX , / , ( 3 ( 2X t E 16 *8 1) 1 

A 

233 

30 

FORMAT (//,15X, 5HD0M l=, E 16. 8, 5X, 4HDAM= ,E 16. 8, 5X, 7HDAMMAX= .E16. 8/ / ) 

A 

234 

31 

FORMAT ( //,40X,11HENGINE DATA, // , 3X.4HN00E ,6X, 4HEN El , 12X, 4HENGJ , 12 

A 

235 


lX,3HENJ,i7X,3HENM, 14X, SHPLY, /, ( IX, 14, IX, E16.8, 2X ,E 16. 8,2X , E 16. 8 , , 2 

A 

236 


2X,ei6.8»2X,E16.8> 1 

A 

257 

32 

FORMAT (//,15X,4HRH0=, E16.8 , lOX, 7HRH0BAR= , E16. 8, // , 1 5X, 6HEIREF- ,E1 

A 

23 8 


16.8,8X,6HEMREF=,E16.a,8X,5HBREF=,El6.8//> 

A 

239 

33 

FORMAT (8(2X, E15.8 » I 

A 

240 

34 

FORMAT l//,20X,17HDIMENSIONAL INPUT , / , 20X, 17 ( IH* 1 , // 1 

A 

241 

35 

FORMAT (//,20X,21HN0N-DIMENSI0NAL INPUT,/, 20X, 21 (IH* »,// 1 

A 

242 

36 

FORMAT {//,1X,9HN00AL PT. , lOX, 2HE I , 16X,2HG J, 16X , 2HEM , 16X , 2HCL, , 16X 

A 

243 


1,2HAC> 

A 

244 

37 

FORMAT ( 5X,I2,4X, E16 .8 , 2X , E 16. 8, 2X, E 16. 8, 2X , FI 6. 8, 2X,E16.6» 

A 

245 

38 

FORMAT (1X,9HN00AL P T. , lOX ,1HA , 17X ,1HB , 17X ,2HS W , 15X , 4HRGSQ> 

A 

246 

39 

FORMAT (//,35X,13HASPECT RAT10= ,E16.8,5X, 19HMA SS R AT 1 0=1 /KAPPA= ,E 1 

A 

247 


16.8//) 

A 

248 

40 

FORMAT (/ ,5X,7HOBEGIN=,2E20.8,5X,7HLBEGIN=,E20.8/) 

A 

249 

41 

FORMAT (7^4x,6HIPROCD,2X,53HITMAX I STOP ISHAPE ISYM ISTART 

A 

250 


1 ILMAX K0PT,/,8{4X,I4),/) 

A 

251 

42 

FORMAT (/,5X,23HESTIMATES OF NAT. FREQ. , 2( 3X ,E 16. 8 ) , / , ( 31 3X, El 6 .8 ) 

A 

252 


D) 

A 

253 

■ 43 

FORMAT (///,20X,23HN0 ENGINES ON THE WINGS,//) 

A 

254 

44 

FORMAT (//,23X,62HELEMENTS OF MATRIX AA (ISTIFF=1, ONE LAMINA PER 

A 

255 


1 MATRIX COLUMN),// ) 

A 

256 

45 

FORMAT (//,40X,21HELEMENTS OF MATRIX AA, / , lOX ,88H( ISTI FF=2 , ONE EN 

A 

257 


IGINE/PYLON OR FUSELAGE PER MATRIX COLUMN, WHEN N0DE=1 FUSELAGE JMP 

A 

258 

— 

RLlED) .// ) - 

A 

259 

45 - 

FORMAT ( / ,30X,14HFLUTTER SPE E0= ,E 1 6. 8/ , 30X , 14( 1H« ) ) . . 

A 

260 

47 

FORMAT ( //.20X.30HNUMBER OF COMPOSTIE LAYERS = 0,/) __ . 

A 

261 

48 

FORMAT ( //,40X,11HENGINE DATA . // , 3X, 4HN0DE ,6X, 4HENE I , I2X.t4HENG J, j2_ 

A 

262 


1X.3HENJ,17X,8HENM( N-1 ) ,9X,3HPYL,/ , ( IX, 14, IX, El 6. 8, 2X ^£16 ._8j ZX^E 16. 

A 

263 

— 

28,2X, E16.8,2X,E16.8) ) _ _ - - 

A 

26 4 

49 

FORMAT (1X,9HN00AL PT. , 1 OX, 3HBCH, 16X , 2HZT ,15X , 3HWT I ) 

A 

265 

50 

FORMAT (5X, I2,6X. 3( E14.6,4X ) ) _ — 

A 

266 

51 

FORMAri //, 1X,110( IH*) ,//) _ __ - - - 

A 

267 


END . - . 

A 

268- 
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SUBROUTINE EVALl (AA,WTI»NL) 

8 

1 


COMPLEX OM,ALtW,E,OBEGIN 

B 

2 


REAL LAMBDA, LBEGIN 

B 

3 

c 


9 

4 

c 

CUMPUTES >*EIGHT OF WING, MASS AND PREPROCESSES El AND GJ 

8 

5 

c 


8 

6 

c 

T^^THICKNESS OF WING COVERS 

8 

7 

c 

TJ=THICKNESS OF EACH LAHINA 

8 

8 

c 

WT=WE1GHT OF WING 

8 

9 

c 

TH£TA=ANGLE OF FILAMENT ORIENTATION 

8 

10 

c 


8 

11 


COMMON // 1X,JX,IERR0R,K0PT/CUM/CM,LAM80A,RH0, EM(60I,EI(61),GJ(61I 

8 

12 


ltA(60)tB(60)fSW(60]»AR,RGS0( 60),MUfE(3t3)fMi60) ,AL(60} ,NNt SWEEP «IP 

8 

13 


2R0C0t ITMAXtl SYMfCH 90) , AC 1 90) /8LKl/RHQBAR,SPAN, BCH (60 1 ,G/BLK2/0 <3. 

B 

14 


33 )»ZT (60) .EIREF,EMREF,OBEGIN, LBEGIN 

B 

15 


DIMENSION AA(8,4), WTI(60) 

8 

16 


IF CNL.EQ.O) RETURN 

B 

17 


03=1. /3. 

8 

18 


SMU= SPAN/ (MU-1) 

8 

19 


WT=.0 

8 

20 


IF IJX.EQ.2) GO TO 1 

8 

21 


PR INT 7 

8 

22 

1 

DO 6 1=1, MU 

B 

23 


X=CI-I)^SMU/SPAN 

B 

24 


X1=(I-2)*SMU/SPAN 

B 

25 


X12=X-.5*SMU/SPAN 

8 

26 


CHORD =8CH-U) 

B 

27 


CH1=CH0RD 

8 

2 8 


IF (I.GT.l) CH1=8CH(!-1) 

B 

29 


B22=.0 

9 

30 


B33=.0 

B 

31 


ZJM12=ZJM1=ZT(I ) 

B 

32 


T=.0 

B 

33 


IF (I.EQ.l) GO TO 2 

B 

34 


ZJM12=(2T(I)+ZTII-1) )*.5 ^ 

B 

35 

2 

T12=.0 

B 

36 


DO 4 IA=1 ,NL 

a 

37 


TJ=AA(l,IA)*(l.+X*AA ( 2f lA) )«-AA(3tIA)*SIN(AA(AfIA)*X) 

B 

38 


TJ12 = AA(1 ♦lA)*U.+Xia*AA(2f lA) )+AA(3f I A ) *S IN ( A A (A, I A ) *X1 2 ) 

B 

39 


THETA=AA( 5f IA)*(1.<-X*AA(6»IA) )+AA(7f IA)*SIN(AA(8,I A)*X) 

ft 

40 


THETA12=AA(5,IA)*(l.+Xl2*AA(6i lA) )tAA(7, IA)*SIN(AA(8f IA)*X12) 

3 

41 


IF (JX.EQ.2) GO TO 3 

3 

42 
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PRINT 8, X,T, THETA 

8 

43 

3 

CALL TRANS ( THETA , All , A33 ) 

8 

44 


ZJ==ZJM1+TJ 

B 

45 


ZJ12^ZJM12+TJ12 

B 

46 


B£TA=03*(ZJ**3-ZJM1**3» 

B 

47 


BETA12=03*(ZJ12**3-'ZJM12^*3) 

B 

48' 

C 

MULTIPLY BY 2 FOR SYMMETRIC SECTION 

8 

49 


B22^B224'2.=«'Ail*BETA 

8 

50 


B33=833+2*<^A35*BETA12 

B 

51 


ZJM1=ZJ . 

B 

52 


ZJM12-ZJ12 

B 

53 


T12=T12+TJ12 

8 

54 

4 

T=T+TJ 

8 

55 

C 

COMPUTE STIFFNESSES OF WING 

B 

56 


EIU )=B22*CH0Rb/EIREF 

8 

57 


CH12^.5*CCH0RD+CH1I 

8 

58 


GJ( I l=4.^833*CH12/EI REF 

8 

59 


IF (I.EQ.l) GO TO 5 

8 

60 


wf=wt^SMU#(CHORD*T+CHl*TP»*RH0&AR*G*5MU#(WTKIItWTI(I-ll I*. 5 

8 

61 

5 

TP=T 

8 

62 

6 

EM( n = CHORD»RHOBAR*2.*T/EMREF*-wtl( n/G/EMREF 

8 

63 

C 

PRINT 1,WT \ 

8 

64 


RETURN ^ i 

8 

65 

C 

' i 

B 

66 

7 

FORMAT {/f lOXf 1HX,15X, 2HTl,15Xf EHTHETAI, 15X, 2HT2,l5X,tHTHETA2, 

B 

67 

8 

FORMAT (3XrE16.8»4<2X,E16. 8)) 

B 

68 


END 

8 

69- 



SUBROUTINE EVAL2 (AA) 

C 

1 


DIMENSION AA(8,4) 

C 

2 


COMMON /CUM/ DM, LAMBDA, RHO,EM(60),EI(61),Gj(61)fA(60),B(60),SH(60) 

C 

3 


i , AR , R G SQ ( 60 ) , MU , E { 3 , 3 ) » W ( 60 ) ♦ A L ( 6 G ) ♦ N N , S V^E EP , i PROC 0 , I TMA X • I S YM , C L ( 

C 

4 


290) ,AC(90) /PYLON/BPYdO) tCPYdO) , DPY (10 ) , ENJ dO ) t ENM ( 10 ) , ENE I ( 10 ) , 

C 

5 


5ENGJ(10)fPYt dO),LOCdOJ,NENG 

C 

6 


COMPLEX OM,W,AL,E 

C 

7 

C 


. c 

8 

C 

PREPROCESSOR OF ENGINE PROPERTIES 

C 

9 

C 


_£L. 

10 . 


IF (NENG.EQ.O) RETURN 

c 

11 


00 1 IJ=1,NENG 

G 

12 


ENEKIJI^AACltlJ) 

C 

13 


ENGJ(IJ)=AA(2dJI 

C 

14 


ENJ(U)=AA(3»IJ) 

c 

15 


ENMd J)=AA(4f I J ) 

C 

16 


PYLC I J)=AA(5,IJ) 

C 

17 


BPy(IJI=PYL( IJ»**3*ENM(IJ)/(3.*ENEHIJJ)/(AR**3) 

c 

18 


CPYIIJ)=ENJ( IJ)«ENM( IJ ) 

c 

19 

"~1 

DPY( IJ)=CPY( IJ )*PYL( 1J)/(ENGJ( IJI^AR) 

C 

20 


RETURN 

G 

21 


END 

C 

22- 
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SUBROUTINE TRANS ( THET A, All t A33 1 

D 

1 


COMMON // IX, JX, lERROR, I0PT/BLK2/D( 3.3) .ZT(60» , EIR EF , EMRE F ,OBEG IN, 

D 

2 


ILBEGIN 

Q 

3 


DIMENSION A{3,3), T(3,3), TI(3,3), ST(3,3) 

D 

4 


CT=COS(THETA) 

0 

5 


SIT=S INI THETA) 

D 

6 


T( 1 / 1 )=:CT«CT 

0 

7 


. TU,2) = SIT*SIT 

0 

8 


T(l,3)=2. ♦CT*SIT 

0 

9 


T(2,l )=TI1,2) 

0 

10 


T(2,2)=T( 1,1) 

D 

11 


T(2,3)=:^T(1,3) 

0 

12 


t(3,l)=~CT*SIT 

D 

13 


TI3.2 )=-TI3,l) 

D 

14 


TI3,3)=T( 1,1)-TI 1,2) 

0 

15 


IF (IX.EQ.2) GO TO 1 

D 

16 


PRINT 8 , ((T(I,J),J=1,5),I=1,3) 

0 

17 

1 

DO 2 1=1,3 

0 

18 


DO 2 J=l,3 

D 

19 


All , J )=0. 

0 

26 

z 

ST(I,J)=.0 

D 

21 


00 3 1=1,3 

0 

22 


DO 3 K=i,3 

D 

23 


DO 3 J=l,3 

0 

24 

3 

STII , K)=STI I ,K)+D( I, J)*T(J,K) 

D 

25 

C 


0 

26 

c 

NOTE TI=TC~THETA) 

0 

27 

c 


D 

28 


TI(l,l)=T(l,i) 

D 

29 


TI (l,2)=T(i,2) 

0 

30 


TI (1,3)=-T{1,3) 

D 

31 


Tn2,I)=T(2,U 

0 

32 


TI(2,2)=T(2,2) 

D 

33 


TI (2,3)=-T(2,3) 

D 

34 


tl(3,l)=-‘TC3,U 

0 

3 5 


Tl(3,2)=-T(3,2) 

D 

36 


Ti(3,2)=T(3,3) 

0 

37 


IF (IX.EQ,2) GO TO 4 

D 

38 


PRINT 8 , I(Tl(I,J),J=l,3),I=l,3) 

0 

39 

4 

CONTINUE 

0 

40 


DO 5 L=l,3 

0 

41 


DO 5 K=l,3 

0 

42 


DO 5 1=1,3 

0 

43 

5 

A(L,K )=A(L,K)+TI(L, I )*ST( I,K ) 

D 

44 


A11=A(1,1)»2. 

0 

45 


A33=A(3,3) 

0 

46 

C 


0 

47 


IF (IX.EQ.2) RETURN 

0 

48 


DO 6 1=1,3 

D 

49 

6 

PRINT 7, (A(I,J), J=l,3) 

0 

50 


RETURN 

0 

51 

C 


0 

52 

7 

FORMAT I/,3(3X,E16.8)) 

D 

5 3 

8 

FORMAT (/,3X,(3(E16.8,2X))) 

0 

54 


END 

0 

5 5- 
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SUBROUTINE VARY ( AA , KV 1, KV2i KVT , DEL . WTI ,NL i I STIFF » 

E 

1 

DIMENSION AA(8,4), WTI160), 

KVU4I, KV2(4», DEL(4» 

E 

2 

DO 1 I=1,KVT 


E 

3 

ri^Kvim 


E 

4 

I2=KV2(I) 


E 

5 

1 AA( 1 1 , 12 >=AA( 1 1, I 2 H-DEt ( 1 1 


E 

6 

IF nSTIFF.EQ*!) CALI EVALl 

IAA,WTI,NL) 

€ 

7 

IF (ISTIFF.EQ.2) CALL EVAL2 

(AA) 

E 

8 

RETURN 


E 

9 

END 


£ 

10- 
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SUBROUTINE RLF lOBEGIN, L8EGIN, I ST ART, ILMAX, IM, IMOR DER ) 

F 

1 

c 

SUBSONIC FLUTTER OF El - GJ WING 

F 

2 


COMPLEX E,GM,H0L0,CDET,CS0ET,DCM,W,AL,WWI60I,AAL(60J ,BB0I,BBI3),QW 

F 

3 


iw,0AL,O,08EGIN 

F 

4 


REAL LAMBOA,LBEGIN,LAM 

F 

5 


DIMENSION X(3I, y(3), 0MH200), 0MRI200), LAM(200), BOAI 5) , ODA 151 

F 

6 


1, IMOROERI5) 

F 

7 


COMMON // IX, JX, IERROR,KOPT/CUM/OM,LAM8DA,RHO, EM(60»,EI(61I,GJC61I 

F 

8 


1,AI60»,B(60I,SWL60I ,AR ,RGSQI 601 ,MU,E(3,31 ,H(601 ,ALl60l,NN, SWEEP, IP 

F 

9 


2R0C0 , I TMA X , I SYM,C L ( 90 1 , AC 1901/ ST AT IC/ 1ST , ALO, CM AC, EL , ISHAP E/P YLON/ 

F 

10 


3BPY( 101,CPY( iOl , DPY( 101 ,ENJ( 101 ,ENM( 101 ,ENEI ClOl , ENG J <101 ,PYLI 101 , 

F 

11 


ALOCdCl ,NENG/FREQ/0M1(51 ,DOMl,OAM,CMMAXl,TdAM,DAMMAX 

F 

12 

c 


F 

13 


NAMELIST /0E8UG1/ XtY 

F 

14 

c 


F 

15 

c 

A=0FFSET OF ELASTIC AXIS FROM SEMI CHORO/SEMI CHORD 

F 

16 

c 

Ak=RECUCEO FREQ.=SEMI CHORD*FREQ. /VELOCITY 

F 

17 

c 

AL=POS. LEADING EDGE UP WING TWIST 

F 

18 

c 

6=SEMI CHORO/WING SPAN 

F 

19 

c 

C- THEOOORSEN LAG FUNCTION 

F 

20 

c 

DYNAMIC PRESSURE=ONE HALF THE AIR DENSITY*VEL. SQD. 

F 

21 

c 

EI=8EN0ING stiffness OF WING/REFERENCE El 

F 

22 

c 

EL=SPAN/FINITE DIFF. SPACING 

F 

23 

c 

EM=MASS PER LENGtH OF SPAN/ REFERENCE MASS PER LENGTH OF SPAN 

F 

24 

c 

GJ=TORSIONAL STIFFNESS OF WING/REFERENCE El 

F 

25 

c 

LAMBDA=2*0YNAMIC PRESSURe*SPAN CU8ED*SEMI CHORO/REF. El 

F 

26 

c 

MU=NO. OF FINITE DIFF. NODES 

F 

27 

c 

RGSO=SQ. OF THE RADIUS 06 GYRATION ABOUT EL. AXIS/SEMI CHORD SQD. 

F 

28 

c 

RH0=PI*SEM1 CHORD SQUARED*AIR OENSITY/REF, MASS PER UNIT SPAN 

F 

29 

c 

SW=OFFSET OF C.G. FROM EL. AXIS/SEMI CHORD 

F 

30 

c 

SWEEP=ANGLE OF WING SWEEP IRADIANSl 

F 

31 

c 

W=POS. UPWARD WING DEFLECTiON 

F 

32 

c 


F 

33 

c 

WING MAY HAVE TWO ENGINES 

F 

34 

c 

ENEI=BENDING STIFFNESS OF P YLON/REFERENCE El 

F 

35 

c 

"ENGS=TQRSIONAL STIFFNESS OF PYLONS/REFERENCE El 

F 

36 

c 

ENJ=TORSIONAL INERTIA OF ENGINE/ ( SPAN*SPAN1 

F 

37 

c 

ENM=ENGINE MASS/WING MASS 

F 

38 

c 

PYL=LENGTH QF PYLON/ SEMI CHORD 

F 

39 


ABSCISA(X3,X2,Xl,Y3,Y2,Yll=X3+( (Y1*Y3-Y2*Y31*( X3-X11 *1 X2-X31 1/1 Yl* 

F 

40 


1Y2*IX1-X2 1+YI*Y3*(X3-X1)+Y2*Y3*(X2-X31 1 ^ ^ “ 

F 

41 


rcL=o 

F 

42 


IF (ISTART.NE.2l GO TO 2 

F 

43 


DM=0BEGIN/(EL*EL1 

F 

44 


LAMB0A=L8EGIN/( EL*EL*EL ) 

F 

45 


IF { LAMBDA. GT. 0. 1 GO TO 1 

F 

46 


IERR0R=4 

F 

47 


RETURN 

F 

48 

1 

HOLD=OM 

F 

49 


0AMMAX1=DAMMAX/(EL*EL*EL1 

F 

50 


GO T O 9 

F 

51 

2 

IM0=1 

F 

^ 


IL=IHORDERUMO) 

F 

53 


DAMMAX1=DAMMAX/ ( EL*EL*EL 1 

F 

54 

3 

IF ( IMO.GT.ILMAXI GO TO 4 

F 

55 


IL=IMCRDERUMO) 

F 

56 


0M2=0M1(IU 

F 

57 


CM2=0M2/(EL*€L1 

F 

58 


GO TO 9 

F 

59 
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4 

If ( ILL.EQ.OJ GO TO 8 


F 

60 


6MIN=B0A( 1) 


F 

61 


IM=l 


F 

62 


IF ( ILMAX.EQ.U GO TO 6 


F 

63 


CO 5 1L=2,ILMAX 


F 

64 


IF (SOAIID.LT.BMIN) i'm=il 


F 

65 

5 

IF (BDA(IL).LT.BMIN J BMIN=BDA(ILI 


F 

66 

6 

00=0DA11M)*EL*EL 


F 

67 


AMDA=BDA( rM)«Et*EL*EL 


F 

68 


print 58 t AMDAfOd. IM 


F 

69 


LAM80A=B0A<I«I 


F 

70 


OMEGA=ODA(IM) 


F 

71 


IF i ISHAPE.EQ.2) GO TO 7 


F 

72 


CALL CFI 


F 

73 


GO TO 35 


F 

74 

7 

IF (KCPT.E0.2I ISTART=2 


F 

75 


IERft0R=0 


F 

76 


LBEGIN=AMDA 


F 

77 


08EGIN=00 


F 

78 


RETURN 


F 

79 

8 

PRINT 59 


F 

80 


IERR0R=1 


F 

81 


RETURN 


F 

.82 

9 

INST A 8=0 


F 

83 


PRINT 50 


F 

84 


IC0UNT=0 


F 

85 


K0UNT=0 


F 

86 


NN=3 


F 

87 


IF <IST.EQ.3) GO TO 42 


F 

88 


IF (I START. EQ. 21 GO TO 17 


F 

89 


GM=CMPLX(0M2*.0) 


F 

90 


HOLO=GM 


F 

91 


00M=CMPLX(DQM1»0. ) 


F 

92 


LAMB0A=0. 


F 

93 


00 10 1=1,3 


F 

94 


IF CI.GE.2) GM=DM+OOM 


F 

95 


CALL CFI 


F 

96 


CALL DEUPPS (E43»CDETt3) 


F 

97 


xm=REAL(COET) 


F 

98 


Y(n=R€AL(OM) 


F 

99 


IF nX.EQ.2) GO TO 10 


F 

100 


PRINT OEBUGl 


F 

101 

10 

CONTINUE 


F 

102 

11 

CISSA=ABSCISA(Y(3> ,YI2) , Y ( 1 ) ,X (3 ) , X ( 2 ) tX(D) 


F 

103 


CM=CMPLX(CISSArO. > 


F 

104 


IF (KGUNT.EQ.O) X0=X(1) 


F 

10 5 


K0UNT=K0UNT+1 


F 

106 


IF «ABS(XMn/ABS(X0I.LE.5.E-A) GO TO 15 


F 

107 


IF (KCUNT.GE.ITMAX) GO TO 14 


F 

108 


CALL CFI 


F 

109 


CALL DEUPPS (£,5,CDET»3I 


F 

110 


CSDET=COET 


F 

111 


Xm = REALCCDETI 


F 

112 


Y(1I=REAL<0M) 


F 

113 


IF MX.EQ.2I GO TO 12 


F 

114 


PRINT DEBUGl 


F 

115 
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12 

CONTINUE 


F 

116 


D0M=D0M/4* 


F 

117 


DO 13 I=:2,3 


F 

118 


CM=OM+OOM 


F 

119 


CALL CFI 


F 

120 


CALL DEUPPS tfe,3,CDET,3) ' “ 


F 

121 


XC I »=REAL(COETI 


F 

122 


Y(I»=REAL(OM» 


F 

123 


IF ( IX.EQ.2r GO TO 13 ~~~ “ 


F 

124 


PRINT DEBUG! ^ 


F 

12 5 

13 

CONTINUE 


F 

126 


GO TO 11 


F 

127 

14 

Ci=QM*EL**2 


F 

128 


AMDAHL AM B D A* EL * ♦ 3 


F 

129 


PRINT 51, 0,AMOA,CSOET,KCUNT 


F 

130 


PRINT 52 


F 

131 


IF UPROCC.E0.2) GO TO A3 


F 

132 

15 

0=0M*EL^*2 


F 

133 


AMDA*LAMBDA«EL*^3 


F 

134 


PRINT 51, 0,AMDA,CDET,KCUNT 


F 

135 


IC0UNT=1 


F 

136 


LAM(1)*0. 


F 

137 


OMR(l)=OM 


F 

138 


GMHlIr^.O 


F 

139 

16 

DAM IN=OAM* U -1 . ) ♦*TNST A8 1 


F 

140 


LAMBOA=LAMBOA+OAMIN 


F 

141 


IF (LAMBDA. Gt.DAMMAXl. AND. ISTART.NE. 2) GO TO 28 


F 

142 

17 

LAMBOA=LAMBDA*COS (SWEEP > *COS (SWEEP) 


F 

143 


ICOUNT-ICOUNT+1 


F 

144 


IF (ICOUNT.GE.ISTOP) GO TO 43 


F 

145 


HOLO=GM 


F 

146 


KOUNT^O 


F 

147 

18 

CALL CFI 


F 

148 


CALL CEUPPS (Ef3,CDET,3) 


F 

149 


CSDET=COET 


F 

150 


AO=REALLCCET) 


F 

151 


AOO-AO 


F 

152 


60=AIMAG(CDET) 


F 

153 


IF { IX^EQo2) GO TO 19 


F 

154 


PRINT 49, COET 


F 

155 

19 

CONTINUE 


F 

156 


GM=0M+CMPLX(D0M1,0*) 


F 

157 


CALL CFI 


F 

158 


CALL DEUPPS (£,3,CDET,3) 


F 

159 


Al-REAL(COET) 


F 

160 


B1=AIMAG(CCET) 


F 

161 


PRWOMR=(Al-AO»/OQMl 


F 

162 


PI WOMR*(Bl-BO)/DOMl 


F 

163 


0M1M=-D0M1 


F 

164 


CM=CM+CMPLX(eMiM,DOMl ) 


F 

165 


CALL CFI 


F 

166 


CALL DEUPPS (E, 3, COET, 3) 


F 

167 


A2=REAL(C0ET> 


F 

168 


B2=AIMAG(CCET) 


F 

169 


PR WO M I - ( A 2- A 0 ) / DO M 1 


F 

170 


PIW0MI=(B2-B0>/D0M1 


F 

171 
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AJAC0BI=PRW0MR*PIW0MI-PRW0MI *PIW0HR 



F 

172 


D0MR= (Bd*PRWOMI-Ad*P IWGMI 1 /AJACOBI 



F 

173 


D0M1*IA0*PIWCMR-B0*PRW0MR)/AJAC0BI 



F 

174 


0M=H0LD+CMPLX(00MR,0GMI I 



F 

175 


H0LD=CM 



F 

176 


IF (I X,EQ,2> GO TO 20 ' 



F 

177 


PRINT 47, PRwdMR,PIMCMR,PRWOMI,PIWOMI 



F 

178 


PRINT 48, DOMR,DOMI,CM 



F 

179 

20 

CONTINUE 



F 

180 


ROM'REAUOMI 



f 

181 


A0M«A1MAG(GMI 



F 



KOUNt^KOUNT+l 



F 

183 


IF {RCM.eO.0.) GO TO 22 



F 

184 


IF (ACM*EQ4»0*) go to 21 



F 

185 


IF IABSIOCHR/ROM|,LT.5.E-3.ANO.A8SIOOMI/AOM).LT.5.E-31 GO TO 24 

F 

186 


GO TO 23 



F 

187 

21 

IF (ABS(0QMR/R0M),LT.l«e-5l GO TO 24 



F 

188 


GO TO 23 



F 

189 

22 

IF <A8SiD“CMI/A0M).LT.l.E-5l GO TO 24 



F 

190 

23 

CONTI NUE 



F 

191 


IF IK0UNT.6E.(2»ITMAXn PRINT 52 



F 

192 


IF IKOUNT.GE. 12*ITMAXI . AND. I PROCO. EQ. 2 1 

GO 

TO 29 

F 

193 


IF (KCUNT.GE.(2*ITMAX) .AND.IPROCD.EQ. 1) 

GO 

TO 24 

F 

194 


GO TO 18 



F 

195 

C 

ROOT LOCUS POINT FOUND 



F 

196 

24 

0=0M*EL**2 



F 

197 


LAMBDA=LAM80A/ (COS (SWEEP )*COS( SWEEP II 



F 

198 


LAMI I CQUNTI -LAMBDA 



F 

199 


CMR(ICOUNT)=ROM 



F 

200 


OMinCOUNT)=AOM 



F 

201 


AMDA=LAMBDA4£L**3 



F 

202 


Print bi, o,amda,c5det,kcunt “ 



F 

203 


IF ( AOM.LT.O..ANd. IC(iUNT.EQ.l) INSTAB*! 



F 

204 


IF (INSTAB*EQ*1> GO TO 27 



F 

205 


IF (K0PT.EQ.2) GO TO 25 



F 

206 

t 




F 

207 

C 

TRACE OUT V - G DIAGRAM 00 NOT STOP AT 1 

FLUTTER 

F 

208 

C 




F 

209 


GO TO 26 



F 

210 

25 

IF CAOM,LT.O..ANO.INSTA8.EQ.O> go to 30 



F 

211 

26 

continCJI 



F 

212 


■MiroilKlfWTWMMiMCTWMM— ■ 



F 

213 


ROM* 2 .♦OMR ( I COUNT I--OMR ( ICOUNT-l i 



F 

214 


A0M=2.*0MI ( ICOUNT I--OMI ( lCOUNT-1 1 



F 

215 


GM*CMPLX(RCM,AOM) 



F 

216 


IF ( rSTART.NE.2I GO TO 16 



F 

217 


IF (I CGUNT.lt. 3) GO TO 16 



F 

218 


CURV=OMI( IC0UNT »-2.*CMI(1C0UNT-1I«-QMI(IC0UNT-2» 

F 

219 


IF (CURV.GT..0.ANC.INSTAB.EQ.0 ) GO TO 44 


F 




IF UC0UNT.GT.12.AN0.ISTART.EQ.2) GO TO 

46 


F 

22 L_. 


GO TO 16 



F 

222 

27 

IF IAOM.GT.O.) go TO 30 



F 

223 


IF (ICOUNT.EQ.l) GO TO 16 



F 

_224 


R0M=2 .*OMR( ICOUNT l-OMRI ICOUNT-l l 



\JL 

??5 


A0M=2.*0MI (ICOUNT »-OMI ( lCOUNT-1 1 



F 

226 


OM*CMPLX(ROM,AOMi 



F 

227 


GO TO 16 



F 

228 
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28 

0AMMAX2=0AMMAX1*EL*EL*EL 

F 

231 


PRINT 55 f OAMMAX2 

F 

232 


IF (IGPT.EQ.2.AN0.ISTART.EQ.2) GOTO 45 

F 

233 

29 

IF (ISTART.EQ.2> RETURN 

F 

234 


BDAM U^OAMMAXl 

F 

235 


ODA( ILI=0. 

F 

236 


IMO=IMO+l 

F 

237 


IF ( ISTART.NE.2I GO TO 3 

F 

238 


GO TO 43 

F 

239 

C 

FLUTTER ENCOUNTERED 

F 

240 

C 

USE ABSCISA TO FIND A BETTER ROOT 

F 

241 

30 

IF ( ICOUNT.lt. 3) GO TO 31 

F 

242 

31 

ROOT=LAM( ICOUNT-1 »-DAMIN*OMI( ICOUNT-1 ) / (OMI ( ICOUNT »-OMI ( ICOUNT- 1) 1 

F 

243 


LAMBOA^ROOT 

F 

244 


IF (ICCUNT.lt. 3) GO TO 33 

F 

245 


00 32 IV=1,3 

F 

246 

32 

Y(IV)=0MR(ICCUNT-3+IV) 

F 

247 


R00T=ABSCISA(Y(3» fY(2) .Ym ,X( 3» ,X(2» fX(in 

F 

248 


GO TO 34 

F 

249 

33 

ROOT*OMRUCOUNT-1 »- ( CMR 1 ICOUNTl -OMR ( ICOUNT-U > ♦OMI ( I COUNT-U/ ( OMI I 

F 

250 


1 ICOUNT) -OMI (ICOUNT-1 )) 

F 

251 

34 

ROMsROOT 

F 

252 


RO=ROI^«EL^EL 

F 

253 


AMDA=LAMBOA*EL*EL*EL 

F 

254 


PRINT 54, AMDAfRO 

F 

255 


IERR0R==0 

F 

256 


LBEGIN=AMOA 

F 

257 


OBEG1N=CMPLX(RO,0.) 

F 

25 8 


IF (I START*EQ^2.AN0.ISHAPE.E0.2) GO TO 42 

F 

259 


IF (1START.EQ.2*AN0.ISHAPE.EQ.1) GO TO 35 

F 

260 


BOACI D ^LAMBOA 

F 

261 


GOAI IL)=RCN 

F 

26 2 


ILL=L 

F 

263 


IM0=IM0>1 

F 

264 


IF II START*NE.2) DAMMAX1=LAM80A+OAM 

F 

265 


IF (ISTART.NE,2,AND*IST.E0.1J GO TO 3 

F 

266 


IF 1 1 SHAPE. EQ. 2) GO TO 42 

F 

267 

C 


F 

268 

C 

MODE SHAPE 

F 

269 

35 

BB0I=l./I E(2,l)*EI3,3)^Er2,3)*E(3,l) ) 

F 

270 


BB(1)=(-E(2,2)*E( 3,3 )+E (3,2 )*E( 2,3)) *3801 

F 

271 


BB(2)=1. 

F 

272 


88 C 3 ) =(-E ( 2 ,1) *E ( 3,2) +E (3 , 1 ) *E ( 2 , 2 ) ) *8801 

F 

273 


LAMBOA^LAMBDA*COS (SWEEP) *00$ (SWEEP) 

F 

274 


00 36 1=1, MU 

F 

275 


J=I 

F 

276 


WW(J)=0. 

F 

277 

36 

AAL( J )=0. 

F 

278 


00 37 NN*1,3 

F 

279 


CALL CFI 

F 

280 


00 37 1=1, MU 

F 

281 


J=I 

F 

282 


WW( J) =WW( JH^BBINN) *W(I ) 

F 

283 
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37 

ML(J)=AAL<J >+B8INN»*AL< n 


F 

284 


PRINT 56 


F 

285 


MUJ=(MU+l»/2 


F 

286 


00 41 


F 

287 


I*J 


F 

288 


IF CREAL(WWIJ)).E0.0..AND.AIMAGCWW(JI),E0.0.) GO TO 38 

F 

289 


CWW=CABS(WW(JI ) 


F 

290 


QWM=CLOG(MUI Ji ) ~ 


F 

291 


PWU-AIMAGICiWWi “ 


F 

292 


GU^O 3^ 


f 

293 

"W 

cww=o. 


F 

294 




F 

295 

39 

IF (REAL<AAL(jn.EQ.0..ANO.AIMAG{AAL(J)) 

.EQ.O. 1 GO TO 40 

F 

296 


CAL»CABS(AAL(jn ~ 


F 

297 


OAL»CLOG(AAL (J ) ) “ 


F 

298 


PAt*AIMAG(QAU 


F 

299 


GO TO 41 


F 

300 

40 

CAL*0. 


F 

301 


PAL=0. 


F 

302 

41 

PRINT 57 f I,CWW,PWW,CAL#PAL 


F 

303 

42 

IF UST.LT,3.ANO.ISTART.E0.2I RETURN 


F 

304 


IF ( IST.tT-3*AN0.ISTART.E0.3) GO TO 7 


F 

30 5 


IF IILL.eQ.2» RETURN 


F 

306 

C 



F 

307 

C 

DIVERGENCE TEST IS PERFORMED ONLY IF IF0R0=2 OR 3 

F 

308 

C 



F 

309 


CALL DTV 


F 

310 


ILL = 2 


F 

311 


IF Cl SHAPE. EQ.l) GO TO 35 


F 

312 


RETURN 


F 

313 

43 

PRINT 60 


F 

314 

C 



F 

315 

44 

IF lOMIMCOUNTI.LT.OMHICOUNT-in GO TO 

16 

F 

316 . 


AM0A=LAM80A«=EL«EL*EL 


F 

317 


PRINT 61, AMOA 


F 

318 

G 



F 

319 

C 

TROUBLE ENCOUNTERED ON A BRANCH CURVE 

CHECK OUT ALL BRANCHES 

F 

320 

C 



F 

321 

45 

ISTART==3 


F 

322 


IERRDR=1 


F 

323 


RETURN 


F 

32 4 

46 

IF (OMKICOUNTJ.LT.OMKICOUNT-in GO TO 

16 

F 

325 


AMOA=LAMBDA^ EL* EL EL 


F 

326 


PRINT 62, AMOA 


F 

327 


IERR0R=2 


F 

328 


ISTART=3 


F 

329 


RETURN 


F 

330 
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47 

FORMAT (/,2X,19HPARTIAL DER IVAT I VE S, / t 4( IX, E16 .8 } » 

F 

333 

48 

FORMAT (3X,5H00MR = ,E16.8,3X,5H00MI=,E16.8,3X,3HriM=,2IE16.8,lXn 

F 

334 

49 

FORMAT </,5X,5HCDET=,2(E16.8,lX) ) 

F 

33 5 

50 

FORMAT (///16X,5H0MEGA23X,6HLAMB0A26X,11H0ETERMINANT23X,13HN0. ITER 

F 

336 


1ATI0NS/5X,9HREAL PARTllX, 14HIMAGINARY PART23X, 9HREAL PART11X,14HI M 

F 

337 


2AGINARY PART/) 

F 

338 

51 

FORMAT (5E20.8,15XI8> 

F 

339 

52 

FORMAT (1X,102HN0 CONVERGENCE AT THE FOLLOWING BRANCH PT. - PROGRA 

F 

340 


IM WILL PROCEED TO NEXT BRANCH PT, PROVIDED IPR0CD=1) 

F 

341 

53 

FORMAT I/, 5X,20HTHE MAGNITUDE OF 0M= , E16.8 , 2H+ I ,E16. 8 ,2X , 8HEXC EEDS 

F 

342 


1 ,E16.8/,5X,37HPR0GRAM WILL PROCEED TO THE NEXT CASE) 

F 

343 

54 

FORMAT C/,30X,42HFLUTTER INSTABILITY ENCOUNTERED AT LAMBOA= , F8 , 5 , 5 

F 

344 


1X» 6H0MEGA= f F8*4/f 30X 

F 

345 


2f/ \ 

F 

346 

55 

FORMAT {/,5X,14HLAMB0A EXCEEDS , FIO* 5 ) 

F 

347 

_5j^ 

FORMAT (//,50X, 17HDEF0RMATI ON SHAP E , // , IX, 8HP0SI TI ON , 25X , 1 OHDEF LEC 

F 

348 


1TION,40X,5HTWI ST,/,23X,9HMAGNlTU0E,20X,aHPHASE,15X,9HMAGNITU0E, 15X 

F 

349 


2,5HPHASE,/) 

F 

350 


FORMAT (3X,I3,10X,E16.8,8X,E16.8,10X,E16.8,8X,E16.8) 

F 

351 

58 

FORMAT (//,10X,62HMINIMUM flutter SPEED FOR PRELIMINARY DESIGN OCC 

F 

352 


lURS AT LAMBDA=,E16,8,10HAND OMEGA=, E16.8, IX, 13H0N BRANCH NO., 12/, 1 

F 

353 


20X,51(1H*»,/) 

F 

354 

59 

FORMAT (//,5X,75HALL BRANCHES SEARCHED OUT BUT NO FLUTTER SPEED FO 

F 

355 


.lUND FOR PRELIMINARY DESIGN) 

F 

356 

“ 60'“ 

FORMAT (//,5Xa6HST0PPE0 ON ISTOP) 

F 

357 

61 

FORMAT (//,5X,73HCURV IS POSITIVE AND INSTAB=0 THEREFORE IEPR0R=1 

F 

358 


1 TEST OUT ALL BRANCHES/ ,5X ,21H0CCURANCE AT L AMBDAN= , E16.8/ ) 

F 

359 


FORMAT (//,5X,74HIC0UNT EXCEEDS 12 WHILE ISTART=2' SET IERR0R=2 AN 

F 

360 


ID SEARCH OUT ALL BRANCHES ,/,5X,20H0CCURANCE AT LAMBDA*, £16.8/) 

F 

361 


END 

F 

362-* 
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SUBROUTINE Cfl 

G 

1 

c 


G 

2 

c 

TO DETERMINE THE REDUCED FLUTTER DETERMINANT 

G 

5 

c 


_JL_ 

4 

c 

ISYM=1, CLAMPED ROOT 

G 

5 

c 

ISYM=2t SYMMETRIC WING FLUTTER AND/OR DIVERGENCE 

G 

6 

c 

ISYM=3f ANTI-SYMMETRIC wING FLUTTER AND/OR DIVERGENCE 

G 

7 


REAL LAMBDA 

G 

8 


COMPLEX CM,W,ALt EtCMZtAK, I If Ct01,D2f QMKlOlf QHdOJ tOALIlOJ ,QMA( 101 

G 

9 


IfOMBI 10» t AQHfAOALf AQALlf0MAl(10» t OMMUI 10» , SH , 5 AL t SHP t AOHl t AK2 1 S ,T , 

G 

10 


2QMA2(l0»,QHia0),QALl(10»,QM2(10»,QMAA,0MBBf01f02f A0,Al,A2tA3,AHl, 

G 

11 


3AH2,AH3fAH4fAALl,AAL2fAAL3,AAL4,W0 

G 

12 


COMMON // IX, JXf IERRCP,KOPT/CUM/OM,LAMBDA,RHOf EM(60» , E II 61 » ,GJ ( 61 » 

G 

13 


l,A(60),B(60)fSWI60)fARrRGSOl60>fMU,EI3f3),WI60l , AL (60) .NN, SWEE P , I P 

G 

14 


2R0CD, ITMAX,ISYM,CL(90), AC(90)/STATlC/IST,AL0,CMACf ELfISHAPE/PYLON/ 

G 

15 


3BPY(10),CPY( 10) ,OPY< 10) ,ENJ(10) ,ENM(10) ,£NE I (1 0 ) ,ENG J ( 10 ) , PYL ( 10 ) , 

G 

16 


4L0C(10) fNENG 

G 

17 


NAMELIST /0E8UG2/ E 

G 

18 


AL0S=CMACS=0, 

G 

19 


IF (I$T«EQ.4I ALOS^ALO 

G 

20 


IF nST*EQ.4) CMACS^CMAC 

G 

21 


BRHO=RHO 

G 

22 


BLAM=LAMBDA 

G 

23 


PI =3. 141592653589793238 

G 

24 


II=CMPCXCO.fl.) 

G 

2 5 


IF nST.GE*2) 0M=0M2 = 0. 

G 

26 


CM2=0M*CM 

G 

2 7 


IF (LAMBDA.EQ. 0, ) GO TO 1 

G 

28 


AK^CM*SQRT TRHO/ ( P I^LAMBDA*AR ) ) 

G 

29 


GO TO 2 

G 

30 

1 

AK=0. 

G 

31 

2 

DO 36 N=1,NN 

G 

32 


MN^O 

G 

33 


AL(l]=ALr21=wm=W(2I^W(3)=0, 

G 

34 


IF C1SYM.EQ.2) GO TO 3 

G 

35 


IF (N.NE.3.AN0.IST.NE.4) W(N<-1)=1. 

G 

36 


IF (N.EQ.3.AN0.IST.NE.4) AL(2)=1. 

G 

37 


IMU=2 

G 

38 


GO TO 4 

G 

39 

3 

IF (N.NE.3) WIN)-U 

G 

40 


IF (N.EQ.3) ALUI-i. 

G 

41 


IMU=1 

G 

42 

4 

CONTINUE 

G 

43 


MU1=MU-1 

G 

44 


MU2=MU-2 _ . . 

G 

45 


MU 3= MU- 3 

G 

46 


CS=COS( SWEEP) 

G 

47 


SN^SINCSWEEP) 

G 

48 


TN=SN/CS 

G 

49 


08=l./8. 

G 

50 


PI 2=2. ♦PI 

G 

51 


AR2=AR*AR 

G 

52 


BJ2=B(1)*B(1 ) 

G 

53 


IF IISYM.NE.3) GO TO 8 

G 

54 

C 

calculate wo 

G 

55 


A0=.5#EI( 1)<-.25*EI(1)*0BAR*CS«-.25*GJ(2)*TN*TN-08*EM(1)*0M2*RGSQ(1I 

G 

56 


1*TN*TN/(AR*AR)+08*EM(1)*SWU )*DBAR*0M2*TN/CS/AR-.5*EI(1)*0BAR2/CS2 

G 

57 


2*.5*E I(1)*DBAR/CS*-.25*EI (2)*DBAR2/CS2-O8*EM(l)*CM2*0BAR2/CS2*O8*EM 

G 

58 


3(1)*SW(1 )*0M2*0BAR/AR/CS*TN 

G 

59 


A2=.5*EI( 1)-.5*EI( l)*0BAR/CS-.25fGJ(2)*TN2*08*EM(l )*0M2*TN2*RGS0(1 

G 

60 


1)/AR2-08*EM(1)*SW(1 )*DBAR*0M2*TN/CS/AR+.5*EI( 1 )*DBAR/CS-. 25*EI (2)* 

G 

61 


2DBAR2/CSfrEI(2)*DBAR/CS-.25*EI(2)*08AR2/CS2>08*EM(l )*QM2*DBAR2/CS2- 

G 

6 2 


308*EM(1)*0M2*SW(1 )*OBAR*TN/AR/CS 

G 

63 


43=-. 5*E 1(2) ♦08AR /C S 

G 

64 
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c 

ADO ON AERODYNAMIC LOADS 

G 

65 


AH1=B{1)*CBC 1I*0M2*RH0-II*CH1I*AR*C*LAMBDA*AKI*0BAB/CS 

G 

66 


AH2 = B Cl ) ♦ B ( 1 ) * ( OM 2* R HO + 1 1 *P I *L AM B 0 A* AK*AR* ( 1 .+C L ( 1 } / P I * { C L C 1 ) / P 1 2* 

G 

67 


lACMA/B(lH*C)+CLm*LAMBDA*AR*B(l»*CI*TN/AR 

G 

68 


AM3=-(CL(l)*C+II*PI*AK*B( in*B(U*LAMBDA*TN 

G 

69 


AHA={CL(1»*(CL( l» /P12*ACMA/B(lH*C-Il*Pl*AK*A(ll »*BJ2*LAMB0A 

G 

70 


AAL1=BJ2/AR*(A(1}*RHC/AR*0M2 + II*PI*AK*LAMB0A*ACMA/PI/BC1I*CU1I*C» 

G 

71 


l*DBAR/C$ ^ ^ ^ 

G 

72 


AAL2=BJ2/AR*(RH0*BJ2/AR*(08+A( 1 l*A Cl 1 /BJ2 l*0M2-2.*PI *I I *AK*LAM80A* 

G 

73 


IBC U*C .5+ACMA*CLC l)*C/PI2/BCin*CCLCl > /PI2<-ACMA/BC 1» I-ACMA/0C1 » *CL 

G 

74 


2(1 )*C*LAMBDA)^TN/AR 

G 

75 


AAL3=-C II *P I*AK*A C 11 -ACMA/B C U*CL Cl » *C» *BJ2/ AR2*LAMBDA*TN 

G 

76 


AAL4=-BJ2*TN/AR2*CII*AK*C08+AC 1 1*A Cl 1/BJ 2 1 *BJ2+AC 1 1 +C I.hACMA/PI^BC 

G 

77 


11)*CL CllfCl^CCLCll/PIB^ACMA/BCin^BCll »*PI*LAMBOA 

G 

78 


A0=A0 + .25*DBAR/CS*C*.5*AH1-.5*AH2-.5*AH3h.5*AH4*TN/ARM-.5*TN/AR*C- 

G 

79 


1.54AAL1-. 5+AAL2-. 5*AAL3+.5*AAL4*TN/AR1 

G 

80 


A2=A2+.25*DBAR/CS*C .5*AH1+.54AH2+.5*AH3-.5*AH4*TN/AR )+. 5*TN/AR» C . 5 

G 

81 


1*AAL1+,5*AAL2+.5*AAL3*.5*AAL4*TN/AR1 

G 

82 


A4=.25*DBAR/CS*AH4+.5<‘TN/AR*AAL4+.5*GJC2)*TN*AR 

G 

83 

t 

AOb ON tNGlN^ t^RM$ ~ ^ 

G 

84 


IF CNENG.EQ.Ol GO TO 7 ~ ^ 

G 

85 


DO 6 IJ=1,NENG 

G 

86 


IF CLCCdJI.NE.l) GO TO 5 

G 

87 


AQAL^CPYI IJ)«0M2«TN/(i.-0M2*DPY( Un 

G 

88 


AQAL1=^CPY( IJ )*0M2/(1 .»0P2*0PY( IJ n 

G 

89 


A0=A0Hh,5*TN/AR*(-AQALM'.25*(-AQALl) 

G 

90 


A2=A2+.5*TN/AR*AQAL+.254AQAL1 

G 

91 

5 

If (L0C(IJI.NE.2) GO TO 6 

G 

92 


AQH=CPY C I J 1 *GM2* ARI-C S/ C 1 .-DPYC IJ 1 1 

G 

93 


AQAL = -'0M2*ENM(IJ) *PYL( IJ) /AR^SN/I l,-BPY(Ijn 

G 

94 


AQAL1=,25*0BAR/AR 

G 

95 


A0=A0+^.25*0BAR/CS*C AQH*AQALH-AQAL*C-AQAL1*TN1*PYL( IJl ) 


96 


A2=A2 + .25*DBAR/CS*C ACH*C-AQAL1 l+AQAL*C AQALl*TN )*PYLC IJ 1 + 1. ) 

G 

97 


A3=A3 + .25*DBAR/CS*C .5*CS/AR«-AQAL*C-.5/AR*SN*PYLC IJ 1 1 1 

G 

98 


A4= A 2 5 *DB AR /C ( A QH^SN ^-AQ AL *CS ) 

G 

99 

6 

CONTINUE 

G 

100 

7 

W0=^(A2*W(2I+A(3)^W(3)+A(4)*AU2n/A0 

G 

101 


wm^OBAR^.S^I W( 2)-W0)/CS 

G 

102 

C 

COMPUTE AERODYNAMIC LOADS AT J==l 

G 

103 


01=C AH1«-AH2+AH3-AH4*TN»*.5*C WC21-W0) +AH4*ALC2I 

G 

104 


Q2=CAAL1+AAL2+AAL3-AAL4*TN)*.5*CWC2)-W01+AAL4*ALC2I 

G 

105 

8 

DO 35 J=IMU,MU 

_G_ 

106 

C 


G 

107 


C=Cl.tl0.61*II*AK*BC J) l/C 1. + 13. 51* II*AK*BCJI1* Cl. +1.774*1 I*AK*BCJl 

G 

108 


11/ C1.+2.745*II*AK*BC Jl) 

G 

109 

c 

PRE CALCULATIONS FOR ENGINE LOADS 

__Q_ 

110 

c 


G 

111 

c 

DETERMINE SLOPES 

G 

112 


IF (J.EOd) GO TO 10 

G 

113 


IF (J.EQ.MU) GO TO 11 

G 

114 


SH=C WCJ+H-WC J-1) )*.5 

G 

115 


SAL--*AL( J-l>*^5 

G 

116 


IF (J*EQ*2) GO TO 9 

G 

117 


SHP=«5*(W(J)-W(J-2)) 

G 

118 


GO TO 12 

G 

119 

9 

SHP=-ALm*AR^TN 


120 


IF (ISYM.EQ.3) SHP=AR«AL( ir/TN 

G 

121 


GO TO 12 

G 

122 

10 

SH=^2.*AR*ALCJ)^TN 

G 

123 


SAL==-AL(1) 

G 

124 


GO TO 12 

G 

125 

11 

SHP = . 5*( WCMUI-W(MUi n 

G 

126 


SH=W(MU)-W(MUU 

G 

127 


SAL=AL(MU)-AL(MU1) 

G 

128 
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12 

CONTINUE 

G 

129 


IF (NENG.EQ.OI GO TO 18 

G 

130 


bo 17 IJ^lfNENG 

G 

131 


IF {J.GT.(,LOC(IJ>+l).OR.J,LT.XLOC(IJI-l)l GO TO 16 

G 

132 


D1=1.-OM2*0PY(IJ) 

G 

133 


D2==1.-0M2*0PYMJ» 

G 

134 

C 


G 

135 


IF (J.NE,LOC(UH GO TO 13 

G 

136 


IF (J.EQ^MU) QMMUCI J) = ENJ(IJ)»ENMMJ)^0M2*AR«( At(jr«SN+SH^CS/AR)*G 

G 

137 


lS/D2-0M2*ENM(IJ»*PYL(rj»/AR*((AL( JI*CS-SH*SN/AR»*PYL ( IJ) i-WI J) ) *SN/ 

G 

138 


201 

G 

139 


QH( IJ » = 0M2<‘ENM( IJ »*( (ALIJ>*CS-SH*SN/AR)*PYL( IJH-W{ jn/Ol 

G 

140 


QAU IJ)=0M2*ENMn J)/(AR*AR)*((AL(J»*CS-SH*SN/AR|*PYL ( IJ » +W( J » » ♦PYL 

G 

141 


KIJ »*CS/01+ENJ(IJ»*ENM( IJ>*0M2*(AU J)*SN*SH*CS/AR)*SN/02 

G 

142 


IF (J.EQ.l) QMK IJ»=-CM2*ENM(IJ)*PYL( IJI/AR*(ALm/CS*PYL( IJ»>H(1» 

G 

143 


1)*SN/D1 

G 

144 


CMA(IJ»-QMB{IJ»=QMAinJ)-QMA2( IJ »*QHl ( IJ »=QAL1 ( IJ»=0M2(IJ 1 = 0. 

G 

145 


GO TO 17 

G 

146 

13 

IF ( J.EQ* aOCnj)-l) ) GO TO 15 

G 

147 


0M8(IJI = ENMC IJr*ENJ( IJ >*0M2*AR«(AU J-ll*SNfSHP^CS/ARl*CS/02-0M2^EN 

G 

148 


1M( IJ ) *PYL ( I J I/ARY ( ( ALU-U*CS-SHP*SN/AR) YPYL ( IJ >+W ( J-1 ) >*SN/6l 

G 

149 


IF (ISYM.NE.3.0R.J.GT.2> GO TO 14 

G 

150 


CHKI JI=0M2YENM( IJ)*DBAR/CS*IW(2)-W0»*,5/D1 

G 

151 


OALK IJ»=0M2*ENM( IJ »*DBAR/AR2*. 5Y( W( 2»-W0)*PYL ( I J) /Ol>CPY ( IJ»*QM2* 

G 

15 2 


l(W(2)-W0»*TN/02 

G 

153 


0M2I IJ)=CPYl IJ)*0M2*AR*(AL1 2»*SNYSH/AR*CS»*CS/D2-aM2*ENM( ij**PYU I 

G 

154 


1J)/AR*U AL(2)YCS-SHYSN/AR)*PYU IJ)+W{ 2) »*SN/01 

G 

155 

14 

CONTINUE 

G 

156 


QMAM J)=QH(IJ>=QALUJI=QWA1( IJ|=QMH IJ l=QMA2(IJ )=0. 

G 

157 


GO TO 17 

G 

158 

15 

OMAA^ENMM J)*ENJ C U)*0M2«AR*CS/02 

G 

159 


QMBB=-0M2YENM(I J»*PYL( I J)/ARYPYU1 J» 

G 

160 


CMA(I J)=QPAA*(-W( JJ l/(2.*AR)*CS-QM88*<W(J»/I 2.YAR ) *SN+W< J+-1)) Y SN/D 

G 

161 


11 

G 

162 


QMA2( IJ>=QMAA>QMBB*SN/Dl 

G 

163 


OMAK IJ»=0M2*ENMMJ)#PYL(I J)*PYUI J»/(2. YARYAR>*SN/01 

G 

164 


QHin j)=QALl (1J )=QM2(U J=.0 

G 

165 


GO TO 17 

G 

166 

16 

OH(I J)=QAUI J) = QMAnj| = QP8( IJ)=QMA1( IJ)=QMA2{1J »=QM1{ IJ»=QMMU( I J} = 

0’ 

167 


lOHK IJJ=QAL1( IJ»=QM2(I J>=0. 

G 

168 

17 

CONTINUE 

G 

169 

c 

AERODYNAMIC LOADS 

G 

170 

— 

BJ2=8C 

G 

171 


ACMA=AC(J)-A(J> 

G 

172 


AQH=B( J)Y(0M2*RH0*B( J)-I1*CL( J)*AR*C*LAMBDAYAK)*W( JH-BJ2*(0M2YRHd* 

G 

173 


i'AIJ )/B(J )+I I*PI*L AMBDA*AK*ARY( 1 .yCL ( J 1 /P I *< CU J ) /PI 2 yACMA /BJ J ) ) YC) 

G 

174 


2+C U J » ylAMBDAYARy B I J » *C) YAL ( J 1 

G 

175 


AQAL= BJ2/ARY(AIJI yRHO/ARYOM 2 + M*PI YAK*LAM8DAYACMA/(B( J» YPI (YGL( Jit 

G 

176. .. 


ICI YW( JI + BJ2/ARY{RH0YBJ2/ARY{08+A( J)yA(J> /BJ2 iYCM2-PI2YMYAKYLAMB0A 

G 

177 


2YB( J > Y{ . 5+ACMAYCC ( J ) /(P 1 2YBI J 1) YC 1 Y( CL( J » API 2+ ACMA AB ( J 11 -ACMA/S ( J 1 

G 

178 


3YCU( J)YCYLAMBDA1YALU> 

G 

179 

e 

ACCOUNT FOR SLOPES IN AERODYNAMIC LOADS 

G 

180 

- 

AQH1=(CLC JIY(CL( J) API2+ACMA/B(J11YC-IIYPIYAKYA(J 11 Y6J2YLAMBDAYTN 

G 

181 


AQH=AQH-(CL(J1YC+IIYPI*AKY8( J) }Y8t JIYLAMBOAYSHYTN+AOHIYSAI. 

G 

182 


AQAL1=-BJ2AC ARYAR) yTNYIM*AKY( 08+A(J)YAUI ABJ2 l YBJ 2+A( J) + ( l.+ACMAA 

G 

183 


1 (P 1 YB( J 1 1 YCL (J 1 YC 1 Y« CL ( J 1 AP 1 2+ACMAAB ( J 1 1 YB( J1 1 Y LAM BO AY PI 

G 

184 


A0AL=AQAL-<IIYPIYAKYA( J)-ACMAAB(J1YCL(J1YC1YBJ2A{ARYAR}YLAMB0AYTNY 

G 

18 5 


iSH^-AOALl*SAL 

G 

186 


65 



APPENDIX F 


c 


G 

187 

c 

PROCEED ALONG SPAN CALCULATING ALCJ«-U AND WCJ+21 

G 

188 


jpi=jn 

G 

189 


JP2*J<-2 

G 

190 



G 

191 


JM2=J-2 

G 

192 


T=£I( JPl) 

G 

193 


S=GJCJP1I 

G 

194 


IF (J,EQ*U GO TO 21 

G 

195 


IF ( J .EQ.2I GO TO 24 

G 

196 


IF CJ.EQ.PU1) GO TO 29 

G* 

197 


IF ( J .EQ.MUl GO TO 33 

G 

198 


S=S+.5*AQAU 

G 

199 


AL( JP1>=-GJ(JI*AL(JM1H-IGJIJP1I+6J( JI-OMaX'EMUJORGSOIJI/IAR^ARn^A 

G 

200 


ILCJ»+CM2*SWIJ»*EM(J>/IAR*AR)*M(J)-AQAL 

G 

201 


W(JP2»=-EIIJM1>*WI JM2» + 2.*IE1( J»frEI( JMII »*W(JMll-(EI ( JPl » (J) 

G 

202 


1+EI ( JR1)-EM( J)*0M2»*W(J )+2.*(E MJPU +EIIJ »)*W( JPl>-EMIjr*SW( J1 *0M2 

G 

203 


2*AL(J)+AQH 

G 

204 

c 


G 

205 

c 

APPLY ENGINE LOADS 

G 

206 


IF INENG.EG*0) GO TO 20 

G 

207 


DO 19 IJ=1*NENG 

G 

208 


T=T+.5«QMA1MJ) 

G 

209 


AL(JP1)=AL( JPli-QALdJI 

G 

210 

19 

W IJ P2 1 =W ( J P 2 ) +QH I IJ ) + . 5 ♦ 1 QMB 1 IJ ) -Q M A ( 1 J ) » 

G 

211 

20 

AL( JP1)=ALIJP1)/S 

G 

212 


W( JP2l=IW(JP2)*.5*AQHl*AL{JPin/T 

G 

213 


GO TO 35 

G 

214 


S=S^-AQAL1*.5 

G 

215 


ALC2» = (GJ(2l-.5*EMm»0M2«RGSQIJI/IAR=»AR) l»AL( 1 )+. 5 *OM 2 *EMm * SW(1 

G 

216 


l)*wm/(AR’«=ARJ-*5^AQAL^2.*ei CU*TN*< (W(2)-W(l) )/AR +A L C 1 1 ♦TN) 

G 

217 


W(3l=2,*(Eim4-EI ( 2n*W( 2>-r2.*Eim+EI( 2>-. 5*0M2«EN(1) )*W(1)-. 5«E 

G 

218 


IM( !)♦ SW< 1 > *0M2* AL ( I » +. 5*AQH+2.*ARA=EI ( 1 l♦TN*AL m 

G 

219 


IF (NENG.EQ.O) GO TO 23 

G 

220 


DO 22 IJ^ltNENG 

G 

221 


T=TvQMA1<IJI^.5 

G 

222 


AL(2)=AL( 2)-QAL(l JI**5+QPi(IJ)^*5/AR*TN 

G 

223 

22 

^ ( 3 ) == W ( 3 ) + . 5 *QH ( 1 J > - . 5 ♦ QM A C IJ I 

G 

224 

23 

ALI2I==ALf 2>/S 

G 

225 


N(3I = (Wm-«5*AQHl*AU2) )/T 

G 

226 


GO TO 35 

G 

227 

24 

S=S+*5«AQAL1 

G 

22 8 


AL ( 3 » =-G J ( 2 » ★AL m ♦ ( G J 1 2 » *G J 1 3 > -EM ( 2 > *0M2* RG SQ 12 )/ ( AR*AR 1 1 * AL ( 2 M- E 

G 

229 


1MI2I*CM2*SW(2»/(AR*ABI*W{2)-AQAL 

G 

23 0 


IF C ISYM.NE.3) GO TO 25 

G 

231 


W(4)-.5^EI{l)«CWC2)-2**wm+W0H2#^En2)«(W(3)-2**W(2KW(in+EI(3) 

G 

232 


1*(W(4»-2.*W(3I+W( 2n-EM(2l*OM2*IWC2>-AU2>*SWI2» l♦•(GJI2l*(AL(2>-^W 

G 

233 


2(2)-W0J/AR*.5*TN)*AR-.5*EM(l »*0M2*(R6SQm*TN/AR-SW(2»*DBAR/CS)<'( W 

G 

234 


3<2l-W0»*.5>*.5*TN/ARd-EIIl> *( W( 2I-2.*WI 1 » + W0) +EI 1 2>* ( WI 31 -2.* W (2» 

G 

235 


4+Wll) l-.5*EMU)*I W(1)-SWI1)*AL(UI*0M2M‘.5*DBAR/CS 

G 

236 


W(4r=W(4)-AQH-.25*TN*AR*Q2-.25*DBAR*01/CS 

G 

237 


WI4)=-W(4I 

G 

238 


GO TO 26 

G 

239 

25 

CONTI NUE 

G 

240 


W(4) = 2«*(Eni)+EII2l )*Wlll-l2.^En i)+4.*EI(21 + ElC3I-0M2*EM<2n*WC2 

G 

241 


1)+2.*IEI (2I + EI(3) l*WI3l-0M2*EM(2l*SWC2)*ALC2)-2.*eim*AR*TN*ALIl» 

G 

242 


2+AQH 

G 

243 
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26 

CONTINUE 

G 

244 


IF (NENG.EQ.OI GO TO 28 

G 

245 


DO 27 IJ=1,NENG 

G 

246 


T=T4^.5*QMAinj) 

G 

247 


ALf 3)=ALC 3HQAUU) 

G 

248 


IF (I SYM.EQ*3) W(4) = W(4) + *25*TN«AR*QAL1( I J ) + .2 5*DB AR /C S* ( QHlil J ) *Q 

G 

249 


1M2( IJ II 

G 

250 

27 

W(4)=W(4)+QH{IJ)-',5«QMAn4) 

G 

251 

28 

AL13)=AL(3>/S 

G 

252 


K(4l=lWf4)+e5*AQHl*AL(3))/T 

G 

253 


GO TO 35 

G 

254 

29 

S=S+«5*AQAL1 

G 

255 


AL(MU )=-GJ(MUll*AL(MU2»*(GJ (MUl)>GJ(MU»-0M2*ENUMUl 1/ {AR*AR 1 *RGSQ{ M 

G 

256 


lUl) »*AL(M01)+EM(MUH*0M2/( AR*AR)*SW(MU1)*W(MU1 »-AQ AC. ' 

G 

257 


IF (NENG.EG.O) GO TO 31 

G 

258 


DO 30 U = 1,NENG 

G 

259 

30 

AL{MU)=AL(MU»-QAL(IJ 1 

G 

260 

31 

AL(MUJ=ALIMU»/S 

G 

261 


E{1»N)=EI (l^U2)*(W (MUl )-2.*W(MU2)+W(MU3J )-2.*EI {MUl >*( W(MU>-2.*W(MU 

G 

262 


11 > t-WC MU2) I-0M2*EM(MU1» *( W( MUll -SW( MUU *AL{ MUl » » -AQH- AQHl * AL ( MU J *.5 

G 

263 


IF (NENG,EQ«0) GO TO 35 

G 

264 


DO 32 IJ=lfNENG 

G 

265 

.32 

e(l,NI=E{l,N>-QH(I J)-.5*(0MB(I J)-,5*IQMAnj)*0MAU IJ » *J 2. *W ( MU >-W{ 

G 

266 


1MU1I))I 

_G_ 

267 


GO TO 35 

_iL 

268 

33 

£<2f N ) = EI (MU1)*(WCMU »-2.*W(MUlH-W(MU2> l-.5*EM( MU)*0M2*( W( MU)-AL(MU 

G 

269 


1)*SW( MU) )-. 5*AQH 

G 

270 


E ( 3 , N) =G J ( MU ) * { AL < MU ) - A L ( MU 1 ) ) - . 5 *EM i MU ) ♦ ( A L ( M U ) *R 6SQ 1 MU ) - S W ( MU ) «= W 

G 

271 


HMU) )’»0M2/(AR*AR)-.5*AQAL 

_G_ 

272 


If (NENG.EQ.O) GO TO 35 

G 

273 


DO 34 IJ=1tNENG 

G 

274 


E<2,N)=E(2,NI-.5*(QH(IJ)+QMB(IJ)+QMMU(IJ) ) 

G 

275 

34 

E(3,N >=£{ 3,NI-.5*QAL(rj| 

G 

276 

35 

CONTINUE 

G 

277 


CONTINUE 

G 

278 


RETURN 

G 

279 

C 


G 

280 


END 

G 

281- 
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SUBROUTINE OEUPPS I A,N/DET,MAX) 


H 

1 


COMPLEX A(MAX,N» f SWAP,DEt,CO,Cl 


H 

2 


COMMON // IX, JX, lERRORfKOPT 


H 

3 * 


C0=(0. ,0.) 


H 

4 




M 

5 


0ET=C1 


H 

6 


NN=N-1 


H 

7 

c 

PIVOT SEARCH 


H 

8 


DO 7 1=1, NN 


H 

9 


CAVM=0, 


H 

10 


11=1 


H 

11 


00 1 J=1I,N 


H 

12 


SWAP=AIJ, in 


H 

13 


AR=ReAL(SwAP) 


H 

14 


AT=AIMAG( SWAP r 


H 

15 


CAVA=ABS ( AR)+ABSC AD 


H 

16 


IF {CAVM.GE.CAVA) GO TO 1 


H 

17 


1R0W=J 


H 

18 


CAVM=CAVA 


H 

19 

1 

CONTINUE 


H 

20 


IF CCAVM.EQ.O.) GO TO 8 


H 

21 

c 

ROW INTERCHANGE 


H 

22 


IF ( IROW.EQ.in GO TO 3 


H 

23 


DET=-DET 


H 

24 


00 2 L=II,N 


H 

25 


SWAP = A(IROW,U 


H 

26 


A{IROh,L»=A( II,U 


H 

27 


A(II,U=SWAP 


H 

28 

2 

CONTINUE 


H 

29 

3 

SWAP= 4(11,11) 


H 

3 0 


CET=DET*SWAP 


H 

31 

C 

NORMALIZE PIVOT ROW 


H 

32 


K=in-i 


H 

33 


DO 4 L=K,N 


H 

34 

A 

A( II, L) = A( II, D/SWAP 


H 

35 

C 

ELIMI NATION 


H 

36 


DO 6 L1=K,N 


H 

37 


SWAP=4(Ll,in 


H 

38 


DO 5 L=K,N 


H 

39 


A(Ll fU=4(LlfL)-A( IlfL »*SWAP 


H 

40 

5 

continue ~~ 


H 

41 

6 

CONTINUE 


H 

42 

7 

CONT INUE 


H 

43 


GO TO 9 


H 

44 

8 

DET=CG 


H 

45 


GO TO 10 


H 

46 

9 

CET=OET*A(N*N) 


H 

A7 

10 

RETURN 


H 

. 


END 


H 

49- 


68 



APPENDIX F 



SUBROUTINE DIV 

I 

1 

c 


1 

? 

c 

CUMPUTES OIVERGENCE SPEED 

I 

3 

c 





REAL LAMBDA 

\ 



DIMENSION X{3)', Y(3I 

J - 

6 


complex E.OMfCOETfWfAL 

1 

7 


COMMON // IX, JXiIERRORtKOPT/CUM/OM, LAMBDA, RHO.EM(60I ,EI(61) ,GJ{61I 

I ’ 

8 


1,A( 60) ,B{ 601 ,SWI60) , AR ,RGSQ(60) ,MU,E(3,3 ) ,W(60 ) , AL <60 ) , NN, SWEEP, IP 

I 

9 


2P0CD, ITMAX,ISYM,CL(90),AC(90)/STATIC/ISt,AL0,CMAC,EL,ISHAPE/PYL0N/ 

1 

10 


3BPYI 10) ,CPY(10) ,DPY(10) ,ENJ(10) ,ENM(10) ,ENEI 1 10 ), ENG J ( 10 ) , PYL 1 10) , 



11 


4LOC(10>,NENG 

I 

12 


ABSCI SA( X3,X2,X1,Y3, Y2,Yl)=X3+( (Y1*Y3-Y2*Y3 )*(X3-X1) *(X2-X3) )/ (Yl* 

I 

13 


iY2*(Xl-X2 ) tYl*Y3* ( X3-Xl) +Y2*Y3*(X2-X3) 1 

1 

14 


IF MST^GE^IO) GO TO 7 

\ ~ 



LAMBDA=0. 

I “ 

16 


KOUNT=0 

i 

17 


IST=2 

I 

18 


NN=3 

f 

19 


PRINT 11 

X 

20 


"TCR ^ ^ ^ ^ ^ ^ ^ ^ ^ “ 

I 

21 


0AM1=0AM 

I 

22 


LAMBDA-DAM! 

I 

23 


CALL CFI 

r 

24 


CALL DEUPP5 IE,3,CDET,3) 

I 

25 


XP=REAL(CCET) 

I 

76 

1 

LAMB0A=LAM8DA+DAMI 

I 

27 


KDUNT=K0UNT+1 

I 

28 


CALL CFI 

I 

29 


CALL OEUPPS (E,3,C0ET,3) 

I 

30 


X(1)=REAL(C0ET) 

I 

31 


ALAM=LAMBOA*EL*EL*EL 

I 

32 


PRINt 10, ALAM,COET,KOUNT 

I 

33 


FMFP=A8S( X( 1I~XP I 

I 

34 


If (FMFP.GE.ABSIXMD.ANC.fMFP.GE.ABSlXP)) GO TO 2 

I 

35 


xp^xcu 

I 

36 


IF (KCUNT.GT.lOO) GO TO 6 

I 

37 


IF ILAMBDA.GT.DAMMAX) GO TO 6 

I 

38 


GO TO 1 

I 

39 

2 

Y( 31 -LAMBDA 

I 

40 


Va )=LAMBDA-DAM1 

I 

41 


LAMBDA=LAMB0A-.5*DAM1 

I 

42 


V<2)=LAM8DA 

I 

43 


CALL CFI 

I 

44 


CALL OEUPPS LE,3»COET»3) 

I 

45 


XI3)==X(1) 

I 

46 


X(il=XP 

I 

47 


X(2I = REAUC0ET) 

I 

48 


XO=X(U 

I 

49" 


ALAM=LAMBDA*EL*EL*EL 

I 

50 


PRINT 10, ALAM,CDET,K0UNT 

1 

51 


AM=YC2I 

I 

52 
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3 LAM6DA=ABSClS<UY(3ltY(2),Y(H,X<3> »X(2>tX<lM I 53 

CALL CFI 1 54 

CALL DEUPPS (E>3,COET,3) I 55 

ALAW=LAMBDA»EL*eL*EL 1 56 

PRINT lOt ALAM,CDET,KOUNT I 57 

X(1) = REAL(CDET) I 58 

YdMLAMBDA I 59 

IP (ABS(LAMBDA-AM).LT.1.E-5) GO TO 5 I 60 

AM=LAMB0A I 61 

IP ( K CUNT. GT. 100) GO TO 6 I 62 

IF (LAMBDA. GT.DAMMAX) GO TO 6 I 63 

' LU=2 I 64 

0AMl=0AMl/4. I 65 

K0UNT=K0UNT»1 I 66 

j no 4 T=LUf3 ^ I 67 

LAMBDA=LAMBDA+0AM1 I 68 

CALL CFI I 69 

CALL DEUPPS (E,3,CDET,3> T 70 

X(I ) = REAL(CDET) i 71 

Y(I)=LAMBCA I 72 

4 CONTINUE I 73 

■ GO TO 3 I 74 

5 ■ ALA«=LAMB0A*EL*EL*EL I 75 

PRINT 8, ALAM,C0ET,XCUNT , ‘ ■' i 76 

RETURN I 77 

6 ALAM=LAMBDA*EL*eL*EL I 78 

PRINT 9, KOUNTf ALAMfCDET I 79 

7" RETURN - - I 80 

TC^l ' '' ' ' ' ' I 81 

8 FORMAT (//|20X,21HD1VERGENCE AT LAMBDA* i E16. 8# 5X, 5HCDET= * 2E16. 8 . 5X I 82 

IflBHNO. OP ITERATI0NS=,I4/) I 83 

9 FORMAT {//,20X,20H01VERGENCE NOT FOUNDf5X,6HKOUNT=,I4,5X,7HLAMBDA= .1 84 

T;Fl6.8<5X,5HCDET=,2E16.e/» 185 

1^ FOR MAT (5X,E16.8,10X,E16.8tlX.E16.8»10X,I2) _ I ^ 

U _ _ forma t (8X,6HLAMaOAt30X.4HCDET,15X.10Hl TERATIONSl I 87 

END I 88- 
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SUBROUTINE NATRL ( ILMAXrOSt OE, OOEL I 

J 

1 

c 


J 

2 

c 

COMPUTES NATURAL fREGUENCIES OF AIRCRAFT 

J 

3 

c 

ILMAX=NUM8ER OF FREQUENCIES TO BE FOUND 

J 

4 

c 

GS=STARTING VALUE FOB SEARCH 

J 

5 

c 

OE^FINAL VALUE FOR SEARCH 

J 

6 

c 

COEL=STEP SIZE TO BE USED IN SEARCH 

J 

7 

c 


J 

8 


REAL LAMBDA 

4 

9 


COMPLEX OM,W,AL,E,COEt 

J 

10 



J 

11 



J 

12 



J 

13 



J 

14 


4LaCU0»,NENG/FREQ/0Ml( 5»fDCMl,DAM,CMMAXltTOAM,OAMMAX 

J 

15 


LAMB0A=0* 

4 



NN=3 

4 

17 


ELD=l./{ EL^ELI 

J 

18 


X=OS*ELO 

J 

19 


XO=OOEL=«'ELD 

J 

20 


XMAX=OE*ELO 

J 

21 


IM=:1 

J 

22 


CM=CMPLX(X, .Oi 

J 

23 


PRINT 7f ISYM 

4 

24 


CALL CFI 

J 

2 5 


CALL OEUPPS (E,3,COET,3) 

J 

26 



J 

27 

1 

x=x4-xn 

J 

28 


OK=CMPLX(X, .0) 

n 

29 


IF (X.GT.XMAXJ GO TO 5 

mm 

30 


CALL CFI 

J 

31 


CALL OEUPPS (E,3,C0ET,3I 

J 

32 


F=REAL(CDET) 

J 

33 


IF UX*EQ.2) GO TO 2 

J 

34 


PRINT a, FP,F,X 

4 

35 

2 

CONTINUE 

J 

36 


FMFP=ABS(F-FP> 

J 

37 


IF (fmfp.ge.A8s«f),ano.fmfp.ge,absifp)i go to 3 

J 

38 


FP=F 

J 

39 


GO TO 1 . 

J 

40 

3 

IF UBSIXD».LT.ABS(X»*.5E-2.AND.XD.GT.0.I GO TO 4 

J 

41 


XD=-»5^XO 

J 

42 


FP=F 

J 

43 



GO TO 1 - - 

4 

44 

4 

CM1(IM)=X*EL*EL 

4 

45 


FP = F 

J 

46 


IF (IM.EQ.ILMAXl RETURN 

J 

47 


IM=IM*1 

J 

48 


XD=ODEL*ELD 

J 

49 


GO TO 1 

J 

50 

5^__ 

PRINT 6f IM,0E 

J 

51 


RETURN 

J 

52 

C 


J 

53 

6 

FORMAT I/,2X,25HSEARCHING FOR ROOT NUMBERt I2.27HNO ROOT FOUND OUT 

J 

54 


ITO OMEGA=,E16.8/) 

J 

55 

T 

FORMAT (/f lX,5HISYM=tI4l 

J 

56 

8 __ 

FORMAT C3X,3HFP=,E16.8»4H F=,E16.8*4H X=,E16.3) 

J 

57 


END 

J 

58- 


71 














APPENDIX F 


Input Data for Sample Problem 


^OPTION ISYM=1,K0PT=2, iSHAPE=l, I LMAX=5, lMURi)ER=l,2/3/4,5$ 

$TRACE iJAM=a, IPR0Ca=2, lTrlAX = lQ, ISTOP=20O/DA.iiMAX=4.,DO,.ll=.l$ 

^.INUUT I X=2, JX=2, I Y=2;J 

ijPLAiM SPAji=20.,A=30*-1.0 2,iJ = 3 0*3.,SW=30*.6,RasU=30*2.r)l,BGH = 30*3.,ZT=30x.05, 
U=5.796E9,1.443E8, .0,1.449E8,5. 79iie3/3*.0,4.32E3,o;JEEP=. 0,lJ=30, 
CL=30*6.2831853,AC=30*-.5$ 

^ERGI NE EiiE i = . 0 , EiilGJ= . 0/ ENJ = , 0/ El'll 1= . 0/ PYL= . 0^ LUC = 0, IJEi'lG=0$ 

99Eo I GR E I KEF = 2 . 4 E8, E,9REF = 2 . 5, RHO= . 0023 7/ RHOliAR=5 . 1, I3REF = 3 . , WT I =30*106.554, 
LiJiJ = 3,i<lL=l,D8=0. ,G=3 2,2$ 

OPARAM ISTI FF=1,AA=.0 3 2,7*.0,KVT=1/KV1=5,I<V2=1,l)EL=. 015 70796,MAX=100$ 

^lilATFi^ 00= . 1 , UE = 6 0 , , UU E L= , 3 0 
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Computer Output for Sample Problem 
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TABLE I.- FLUTTER SPEED COMPARISONS FOR UNIFORM 
CANTILEVER RECTANGULAR STRAIGHT WING ’ 

OF ASPECT RATIO 6.67 


Program 

or 

reference 

Description 

Flutter speed, 
km/hr 

Percent 
difference 
with exact 
solution 

Flutter 

frequency, 

Hz 

4 and 5 

Exact analysis 

494 

— 

11.25 

COMBOF 

10 finite-difference stations 

484.2 

-2.0 

11.20 

COMBOF 

15 finite -difference stations 

483.4 

-2.1 

11.24 

COMBOF 

20 finite -difference stations 

483.1 

-2.2 

11.26 

COMBOF 

25 f inite -difference stations 

483.1 

-2.2 

11.27 

SADSAM, ref. 7 

1 finite element 

447.1 

-9.5 


SADSAM, ref. 7 

10 finite elements 

472.5 

-4.4 



TABLE n.- FLUTTER SPEED COMPARISONS FOR THE SYMMETRIC 
FLUTTER OF UNIFORM RECTANGULAR STRAIGHT WING 
OF ASPECT RATIO 6.67 AND WITH ATTACHED 
FUSELAGE AND TIP WEIGHTS 


Program 

or 

reference 

Case 

(a) 

Description 

Flutter speed, 
km/hr 

Flutter 

frequency, 

Hz 

Branch 

5 

1 

Exact analysis 

1054 

3.05 

First torsion 

COMBOF 

1 

31 finite -difference stations 

1049 

3.05 

First torsion 

SADSAM, ref. 7 

1 

10 finite elements 

1049 

3.01 

First torsion 

COMBOF 

1 

31 finite -diff erence stations 

994 

1.81 

First bending 

SADSAM 

1 

10 finite elements 

990 

1.81 

First bending 

COMBOF 

1 

31 finite -difference stations 

943 

16.48 

Second bending 

SADSAM 

1 

10 finite elements 

961 

16.2 

Second bending 

5 

2 

Exact analysis 

826 

3.05 

First torsion 

COMBOF 

2 

31 finite-difference stations 

821 

3.02 

First torsion 

SADSAM 

2 

10 finite elements 

837 

3.00 

First torsion 


^Case 1: center of gravity of tip weights coincides with elastic axis 
Case 2; center of gravity pf tip weights located aft of elastic axis 
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TABLE m.- FLUTTER SPEED COMPARISONS FOR UNIFORM 
CANTILEVER RECTANGULAR SWEPT WING 
OF ASPECT RATIO 12.4 


Sweep angle, 
deg 

Flutter speed, km/hr 

Percent 

difference 

with 

experiment 

Reference 6 

COMBOF 

Experimental 

Analytical 

0 

368 

346 

354 

-3.8 

30 

378 

368 

370 

-2.1 

45 

433 

426 

423 

-2.3 

60 

563 

568 

541 

-2.1 


TABLE IV.- PROPERTIES OF BORON-EPOXY COMPOSITE, UNIFORM, 
CANTILEVERED, RECTANGULAR STRAIGHT WING 
OF ASPECT RATIO 6.67 


Elastic moduli: 

E^ = 276 GN/m2 

E^ = 27.6 GN/m2 

v^t = 0.25 

G^^ =10.3 GN/m2 


Wing properties: 
a = -0.34 

b = 1.0 

= 0.2 

r^2 = 0.29 
a 

C, = 2tt 
’'a 

a^ - -0-5 
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Lumped fuselage 
mass aft 












gw 



Figure 3.- Laminated, balanced ply, filamentary composite wii^ box. 
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Figure 4.- Flow diagram of COMBOF. 
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Angle of filamem orientation, 6 , deg 

Figure 6.- Effect of filament orientation on flutter speed of 
several swept cantilever wings for m^m^K = 8. 
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Flutter parameter. 


Unstable 


✓ / 

f / 

/ 

/ 

/ 

/ 

/ 

/ / 


Stable 


Angle of filament orientation, d, deg 

Figure 7. - Effect of filament orientation on the flutter speed of 
uniform rectangular straight cantilever wing with aspect 
ratio of 6.67 and m/m^K = 64. 
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Flutter parameter, X 

Figure 8.- Projection of root locus branches on imaginary plane at filament orientations near discontinuity 
of lowest flutter speed of uniform rectai^ular straight wing of aspect ratio 6.67 and m/mj./c = 64. 



Flutter parameter 





Flutter parameter, 







1. Report No. 

NASA TN D-7539 


2, Government Accession No. 


4. Title and Subtitle 

FLUTTER ANALYSIS OF SWEPT-WING SUBSONIC AIRCRAFT 
WITH PARAMETER STUDIES OF COMPOSITE WINGS 


7. Author{s) 

J. M. Housner and Manuel Stein 


9. Performing Organization Name and Address 

NASA Langley Research Center 
Hampton, Va. 23665 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 


15, Supplementary Notes 


3. Recipient's Catalog No. 


5. Report Date 

September 1974 


6. Performing Organization Code 


8. Performing Organization Report No. 

L-9260 


10. Work Unit No. 

501-22-03-06 


1 1 . Contract or Grant No. 


13. Type of Report and Period Covered 
Technical Note 


14. Sponsoring Agency Code 


16. Abstract 

This paper presents a computer program which has been developed for the flutter 
analysis, including the effects of rigid-body roll, pitch, and plunge, of swept -wing subsonic 
aircraft with a flexible fuselage and engines mounted on flexible pylons. The program 
utilizes a direct flutter solution in which the flutter determinant is derived by using finite 
differences, and the root locus branches of the determinant are searched for the lowest 
flutter speed. In addition, a preprocessing subroutine is included which evaluates the 
variable bending and twisting stiffness properties of the wing by using a laminated, balanced 
ply, filamentary composite plate theory. The program has been substantiated by compari- 
sons with existing flutter solutions. 

The program has been applied to parameter studies which examine the effect of fila- 
ment orientation upon the flutter behavior of wings belonging to the following three classes: 
wings having different angles of sweep, wings having different mass ratios, and wings 
having variable skin thicknesses. These studies demonstrated that the program can per- 
form a complete parameter study in one computer run. The program is designed to detect 
abrupt changes in the lowest flutter speed and mode shape as the parameters are varied. 


17. Key Words (Suggested by Author(s)) 
Subsonic flutter 
Composite wing 


18. Distribution Statement 

Unclassified -- Unlimited 


19. Security Classif. (of this report) 
Unclassified 


20. Security Classif. (of this page) 

21 . No. of Pages 

Unclassified 

106 


STAR Category 32 


22. Price* 

$4.00 


For sale by the National Technical Information Service, Springfield, Virginia 22151 
























NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
WASHINGTON, D,C. 20546 


OFFICIAL BUSINESS 

PENAL.TY FOR PRIVATE USE $300 SPECIAL FOURTH-CLASS RATE 

BOOK 


POSTAGE AND FEES PAID 
NATIONAL AERONAUTICS AND 
SPACE ADMINISTRATION 
451 



POSTMASTER : 


If Undeliverable (Section 158 
Postal Mannal ) Do Not Return 


aeronautical and space activities of the United States shall be 
conducted so as to contribute to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof/' 

— ^National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientifc «n<l 

technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 

Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. Also includes conf erence 
proceedings with either limited or unlimited 
distribution. 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include final reports of major 
projects, monographs, data compilations, 
handbooks, sourcebooks, and special 
bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may he of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs, 
Technology Utilization Reports and 
Technology Surveys. 


Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION OFFICE 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 




